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Abstract 

This thesis responds to the challenges of using a large number, such as thousands, of features 
in regression and classification problems. 

There are two situations where such high dimensional features arise. One is when high 
dimensional measurements are available, for example, gene expression data produced by 
microarray techniques. For computational or other reasons, people may select only a small 
subset of features when modelling such data, by looking at how relevant the features are to 
predicting the response, based on some measure such as correlation with the response in the 
training data. Although it is used very commonly, this procedure will make the response 
appear more predictable than it actually is. In Chapter 2, we propose a Bayesian method 
to avoid this selection bias, with application to naive Bayes models and mixture models. 

High dimensional features also arise when we consider high-order interactions. The num- 
ber of parameters will increase exponentially with the order considered. In Chapter 3, we 
propose a method for compressing a group of parameters into a single one, by exploiting 
the fact that many predictor variables derived from high-order interactions have the same 
values for all the training cases. The number of compressed parameters may have converged 
before considering the highest possible order. We apply this compression method to logistic 
sequence prediction models and logistic classification models. 

We use both simulated data and real data to test our methods in both chapters. 
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Chapter 1 
Introduction 



1.1 Classification and Regression 

Methods for predicting a response variable y given a set of features x = (x\, . . . ,x p ) are 
needed in numerous scientific and industrial fields. A doctor wants to diagnose whether a 
patient has a certain kind of disease from some laboratory measurements on this patient; 
a post office wants to use a machine to recognize the digits and characters on envelopes; a 
librarian wants to classify documents using a pre-specified list of topics; a businessman wants 
to know how likely a person is to be interested in a new product according to this person's 
expenditure history; people want to know the temperature tomorrow given the meteorologic 
data in the past; etc. Many such problems can be summarized as finding a predictive function 
C linking the features a: to a prediction for y: 



The choice of function C depends also on the choice of loss function one wishes to use in 
making a decision. In scientific discussion, we focus on finding a probabilistic predictive 
distribution: 





1 



1 Introduction 



2 



P(y | x) (1.2) 

Here, P(y \ x) could be either a probability density function for continuous y (a regression 
model), or a probability mass function for discrete or categorical y (a classification model). 
Given a loss function, one can derive the predictive function C from the predictive dis- 
tribution P(y | x) by minimizing the average loss in the future. For example, when y is 
continuous, if we use a squared loss function L(y,y) = (y — y) 2 , the best guess of y is the 
mean of P(y | a;); if we use an absolute loss function L(y, y) = \y — y\, the best guess is the 
median of P(y \ x); and when y is discrete, if we use — 1 loss function L(y, y) = I(y ^ y), 
the best guess is the mode of P(y \ x). 

One approach to finding P(y | x) is to learn from empirical data — data on a num- 
ber of subjects that have known values of the response and values of features, denoted by 
{(y (1) ,a; (1) ), . . . , (y( n \x {n) )}, or collectively by (y train , aj train ). This is often called "training" 
data, and the subjects are called "training" cases, as we are going to use these data to "train" 
an initially "unskilled" predictive model, as discussed later. In contrast, a subject whose re- 
sponse and features are denoted by (y*,x*), for which we need to predict the response, is 
called a "test" case, because we can use the prediction result to test how good a predictive 
model is if we are later given the true y*. 

There are many methods to learn from the training data (Hastie, Tibshirani and Fried- 
man 2001 and Bishop 2006). One may estimate P(y* \ x*) using the empirical distribution of 
the responses in the neighbourhood of x* in some metric, as in the /c-nearest-neighbourhood 
method. Such methods are called nonparametric methods. In this thesis, we consider para- 
metric methods, in which we use a closed-form function with unknown parameters to model 
the data. Once the parameters are inferred from the training data we can discard the train- 
ing data because we only need the parameters of the "trained" model for making predictions 
on test cases. 
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One class of parametric methods, called conditional modelling methods, start by defining 
P(y | a;) as a function involving some unknown parameters, denoted by 6. These parameters 
will be inferred from training data. For continuous y, the simplest and most commonly used 
form for P(y \ x) is a Gaussian model: 

P(y | x, (3, a) = -j= exp I — I (1.3) 

For a discrete y that takes K possible values 0, . . . , K — 1, one may use a logistic form for 
P(y | x): 

P{y = k | *, 6) = (L4) 

The function f(x,(3) or functions fj(x,{3j) link x to y. They are often linear functions of 
x, but may be also nonlinear functions of x defined, for example, by multilayer perceptron 
networks. Our work in Chapter 3 uses linear logistic models. 

Another class of methods model the joint distribution of y and x by some formula with 
unknown parameters 6, written as P(y,x | 6). The conditional probability P(y \ x,6) can 
be found by: 



P(y, x I 0) , , 



Examples of such P(y, x, 0) include naive Bayes models, mixture models, Bayesian networks, 
and Markov random fields, etc., all of which use conditional independency in specifying 
P(y, x | 0). For example naive Bayes models assume all features x are independent given y. 
Our work in Chapter 2 uses naive Bayes models and mixture models. 

There are two generally applicable approaches for inferring 6 from the training data. One 
is to estimate 6 using a single value, 0, that maximizes the likelihood function or a penalized 
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likelihood function, i.e., the value that best fits the training data subject to some constraint. 
This single estimate will be plugged in to P(y | x, 6) or P(y, x \ 0) to obtain the predictive 
distribution P(y* \ x*,6) for a test case. 

Alternatively, we can use a Bayesian approach, in which we first define a prior distribu- 
tion, P(0), for 0, which reflects our "rough" knowledge about 6 before seeing the data, and 
then update our knowledge about 6 after we see the data, still expressed with a probability 
distribution, using Bayes formula: 

P(g| ^ ^^^iw (16) 

P(0 I y train ; ajt^in) j s called the posterior distribution of 6. The joint distribution of a test 
case (y*, x*) given the training data (y train ; cc train ) is found by integrating over 8 with respect 
to the posterior distribution: 

P(y*,x* | y train , a; train ) = y P(y*,x* | y tia[n ,x t ^ in ,0)P(0 \ y train , a; train ) rf0 (1.7) 

The predictive distribution can then be found as P(y*, as* | y' 1 ^ 1 ^ x tI!iin )/P(x* \ y train 7 a3 train ), 
which will be used to make predictions on test cases in conjunction with our loss function. 



1.2 Challenges of Using High Dimensional Features 

In many regression and classification problems, a large number of features are available for 
possible use. DNA microarray techniques can simultaneously measure the expression levels 
of thousands of genes (Alon et.al. 1999, Khan et.al. 2001); the HIRIS instrument for the 
Earth Observing System generates image data in 192 spectral bands simultaneously (Lee 
and Landgrebe et.al. 1993); one may consider numerous high-order interactions of discrete 
features; etc. 
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There are several non-statistical difficulties in using high-dimensional features, depend- 
ing on the purpose the data is used for. The primary one is computation time. Models for 
high dimensional data will require high dimensional parameters. Consequently, the time for 
training the model and making predictions on test cases may be intolerable. For example, a 
speech recognition program or data compression program must be able to give out the pre- 
diction very quickly to be practically useful. Also, in some cases, measuring high dimensional 
features takes substantially more time or money. 

Serious statistical problems also arise with high dimensional features. When the number 
of features is larger than the number of training cases, the usual estimate of the covariance 
matrix of features is singular, and therefore can not be used to compute the density function. 
Regularization methods that shrink the estimation to a diagonal matrix have been proposed 
in the literature (Friedman 1998, Tadjudin and Landgrebe 1998, 1999). Such methods usually 
need to adjust some parameters that control the degree of shrinkage to a diagonal matrix, 
which may be difficult to determine. Another aspect of this problem is that even a simple 
model, such as a linear model, will overfit data with high dimensional features. Linear logistic 
models with the coefficients estimated by the maximum likelihood method will have some 
coefficients equal to oo; the solution is also not unique. This is because the training cases can 
be divided by some hyperplanes in the space of features into groups such that all the cases 
with the same response are in a group; indeed, there are infinitely many such hyperplanes. 
The resulting classification rule works perfectly on the training data but may perform poorly 
on the test data. Overfitting problems usually arises because one uses more complex models 
than the data can support. For example, when one uses a polynomial function of degree n 
to fit the relationship between x and y in n data points (y^\x^), there are infinitely many 
such polynomial functions that go exactly through each of these n points. 

A sophisticated Bayesian method can overcome the overfitting problem by using a prior 
that favours simpler models. But unless one can analytically integrate with respect to the 
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posterior distribution, implementating such a Bayesian method by Markov chain sampling is 
difficult. With more parameters, a Markov chain sampler will take longer for each iteration 
and require more memory. It may need more iterations to converge, or get trapped more 
easily in local modes. Also, with high dimensional features, it is harder to come up with a 
prior that reflects all of our knowledge of the problem. 

1.3 Two Problems Addressed in this Thesis 

For the above reasons, people often use some methods to reduce the dimension of features 
before applying regression or classification methods. However, a simple implementation of 
such a "preprocessing" procedure may be invalid. For example, we may first select a small 
set of features that are most correlated with the response in the training data, then use 
these features to construct a predictive distribution. This procedure will make the response 
variable appear more predictable than it actually is. This overconfidence may be more 
pronounced when there are more features available, as more actually useless features will 
by chance pass the selection process, especially when very few useful features exist. In 
Chapter 2, we propose a method to avoid this problem with feature selection in a Bayesian 
framework. In constructing the posterior distribution of parameters, we condition not only 
on the retained features, but also on the information that a number of features are discarded 
because of their weak correlations with the response. The key point in our solution is that 
we need only calculate the probability that one feature is discarded, then raise it to the 
power of the number of discarded features. We therefore can save much computation time 
by selecting only a very small number of features for use, and at the same time make well- 
calibrated predictions for test cases. We apply this method to naive Bayes models and 
mixture models for binary data. 

A huge number of parameters will arise when we consider very high order interactions of 
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discrete features. But many interaction patterns are expressed by the same training cases. 
In Chapter 3, we use this fact to reduce the number of parameters by effectively compressing 
a group of parameters into a single one. After compressing the parameters, there are many 
fewer parameters involved in the Markov chain sampling. The original parameters can later 
be recovered efficiently by sampling from a splitting distribution. We can therefore consider 
very high order interactions in a reasonable amount of time. We apply this compression 
method to logistic sequence prediction models and logistic classification models. 



1.4 Comments on the Bayesian Approach 

The Bayesian approach is sometimes criticized for its use of prior distributions. Many people 
view the choice of prior as arbitrary because it is subjective. The prior is the distribution 
of that generates, through a defined sampling distribution, the class of data sets that will 
enter our analysis. Thus, there is only one prior that accurately defines the characteristics 
of the class of data sets, which may be described in another way, such as in words. Different 
individuals may define different classes of data sets. The choice of prior is therefore subjec- 
tive, but not arbitrary, since we may indeed decide that a prior distribution is wrong if the 
data sets it generates contradict our beliefs. Typically we choose a diffuse prior to include 
a wide class of data sets, but de-emphasize some data sets we believe less likely to appear 
in our analysis, for example a data set generated by a linear logistic model with coefficient 
equal to 10000 for a binary feature. This distribution is therefore also phrased as expressing 
our prior belief, or our "rough" knowledge about which 6 may have generated our data set. 

There is usually useful prior information available for a problem before seeing any data set, 
such as relationships between the parameters (or data). For example, a set of body features 
of a human should be closer to those of a monkey than to other animals. A sophisticated prior 
distribution can be used to capture such relationships. For example, we can assign the two 
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groups of parameters, which are used to define the distribution of body features of a human 
and a monkey, a joint prior distribution in which they are positively correlated (Gelman, 
Bois and Jiang 1996). We usually construct such joint distributions by introducing some 
extra parameters that are shared by a group of parameters, which may also have meaningful 
interpretations. One way is to define the priors of the parameters of likelihood function in 
terms of some unknown hyperparameters, which is again given a higher level distribution. For 
example, in Automatic Relevance Determination (ARD) priors for neural network regression 
(Neal 1996), all the coefficients related to a feature are controlled by a common standard 
deviation. Such priors enable the models to decide whether a feature is useful automatically, 
through adjusting the posterior distribution of the common standard deviation. Similarly, 
in the priors for the models in Chapter 2 we use a parameter a to control the overall degree 
of relationship between the features and response. Our method for avoiding the bias from 
feature selection has the effect of adjusting the posterior distribution of a to be closer to the 
right one (as would be obtained using the complete data), by conditioning on all information 
known to us, both the retained features and information about the feature selection process. 
Another way of introducing dependency is to express a group of parameters as the functions 
of a group of "brick" parameters. For example, in Chapter 3, the regression coefficients 
for the highest order interaction patterns are expressed as sums of parameters representing 
the effects of lower order interaction patterns. Such priors enable the models to choose the 
orders automatically. 

Once we have assigned an appropriate prior distribution for a problem, all forms of 
inference for unknown quantities, including the unknown parameters, can be carried out very 
straightforwardly in theory using only the rules of probability, since the result of inference 
is also expressed by a probability distribution. These predictions are found by averaging 
over all sets of values of 6 that are plausible in light of the training data. Compared with 
non-Bayesian methods, which use only a single set of parameters, the Bayesian approach has 
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the following advantages from a practical viewpoint. 

First, the prediction is automatically accompanied by information on its uncertainty in 
making predictions, since the prediction is expressed by a probability distribution. 

Second, Bayesian prediction may be better than prediction based on only a single set of 
parameters. If the set of parameters that best explains the training data, such the MLE, 
is not the true set of parameters that generates the training data, we still have the chance 
to make good predictions, since the true set of parameters should be plausible given the 
training data and therefore will be considered as well in Bayesian prediction. 

Third, sophisticated Bayesian models, as described earlier, will self- adjust the complexity 
of a model in light of the data. We can define a model through a diffuse prior that can cover 
a wide class of data sets, from those with a low level of complexity to those with a high 
level of complexity. If the training data does not favour the high complexity, the posterior 
distribution will choose to use the simple model. In theory we do not need to change 
the complexity of a model according to the properties of the data, such as the number of 
observations. The overfitting problem in applying a complex model to a data set of small 
size is therefore overcome in Bayesian framework. Although more complex models may make 
the computation harder, Bayesian methods are, at least, much less sensitive to the choice of 
model complexity level than non-Bayesian methods. 

Bayesian inference, however, is difficult to carry out, primarily for computational reasons. 
The posterior distribution is often on a high dimensional space, often takes a very complicated 
form, and may have a lot of isolated modes. Markov chain Monte Carlo (MCMC) methods 
(Neal 1993, Liu 2001 and the references therein) are so far the only feasible methods to draw 
samples from a posterior distribution (Tierney 1994). In the next section, we will briefly 
introduce these methods. However, for naive Bayes models in Chapter 2 we do not use 
MCMC, due to the simplicity of naive Bayes models. 
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1.5 Markov Chain Monte Carlo Methods 

We can simulate a Markov chain governed by a transition distribution T(0' \ 0) to draw 
samples from a distribution 7r(0), where e S, if T leaves ir invariant: 

/ tt(0)T(0' \0)d0 = ir(0') (1.8) 

and satisfies the following conditions: the Markov chain should be aperiodic, i.e., it does 
not explore the space in a cyclic way, and the Markov chain should be irreducible, i.e., the 
Markov chain can explore the whole space starting from any point. Given these conditions, 
it can be shown there is only one distribution ir satisfying the invariance condition (1.8) for 
a Markov chain transition T if there is one. (The condition of aperiodicity is not actually 
required for Monte Carlo estimation, but it is convenient in practice if a Markov chain is 
aperiodic, since we have more freedom in choosing the iterations for making Monte Carlo 
estimation. And it is obviously required to ensure that the result in (1.9) is true.) 

Let us denote a Markov chain by 0^,0^,... . (Roberts and Rosenthal 2004) shows 
that if a Markov chain transition T satisfies all the above conditions with respect to ir, then 
starting from any point 0q for 0^°\ the distribution of 0^ will converge to ir: 

lim P(0 (n) = | (o) = O ) = tt(0), for any 0,0 o eS (1.9) 

n— >oo 

In words, after we run a Markov chain sufficiently long, the distribution of 0^ (regardless 
the starting point) will be close to the target distribution vr(0) in some metric (Rosenthal 
1995 and the references therein). We can therefore use the states afterward as samples from 
7r(0) (though correlated) for making Monte Carlo estimations. It is tremendously difficult 
to determine in advance how long we should run for an arbitrary Markov chain, though we 
can do this for some types of Markov chains (Rosenthal 1995). In practice we check the 
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convergence by running multiple chains starting from different points and see whether they 
have mixed at a certain time (see for example Cowles and Carlin 1996, and the references 
therein). 

It is usually not difficult to construct a Markov chain transition T that satisfies the 
invariance condition for a desired distribution n and the other two conditions as well, based 
on the following facts. First, one can show that a Markov chain transition T leaves tt invariant 
if it is reversible with respect to ir: 

tt(0)T(0' I 0) =tt(0')T(0 I 0'), for any 0', G S (1.10) 

It therefore suffices to devise a Markov chain that is reversible with respect to tt. Second, 
applying a series of Markov chain transition Tj that have been shown to leave 7r invariant 
will also leave ir invariant. Also, applying a series of appropriate Markov chain transition Tj 
that explores only a subset of S can explore the whole space, S. 

Gibbs sampling method (Geman and Geman 1984, and Gelfand and Smith 1990) and 
the Metropolis-Hastings method (Metropolis et. al. 1953, and Hastings 1970) are two basic 
methods to devise a Markov chain transition that leaves ir invariant. We usually use a 
combination of them to devise a Markov chain transition satisfying the above conditions for 
a complicated target distribution ir. 

Let us write = (#i, . . . , 9 P ). Gibbs sampling defines the transition from ( -* _1 * ) to 6^ as 
follows: 
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Draw 9{ from tt{6\ 
Draw $2 from ^(#2 



At) At-i) At-1)^ 
V\ , t>3 , • • • , J 



Draw 6\ from 7r(^ 



/)(*) At) At-i) At-i)\ 



Draw 9 P from n(9. 



p 




The order of updating 6i can be any permutation of 1, ... ,p. One can show each updating 
of 9i is reversible with respect to ir(0), and a complete updating of all 9% therefore leaves 
7r(0) invariant. Sampling from the conditional distribution for 9i can also be replaced with 
any transition that leaves the conditional distribution invariant, for example, a Metropolis- 
Hastings transition as described next. 

The Metropolis-Hastings method first samples from a proposal distribution T(0* \ O^" 1 ^) 



to propose a candidate 6*, then draws a random number U from the uniform distribution 
over (0,1). If 



we let 0^ ' = 6*, otherwise we let 0^ ' = 6^ ' . One can show that such a transition is 
reversible with respect to it, and hence leave 7r invariant. 

1.6 Outline of the Remainder of the Thesis 

We will discuss in detail our method for avoiding bias from feature selection in Chapter 
2, with application to naive Bayes models and mixture models. In Chapter 3 we discuss 
how to compress the parameters in Bayesian regression and classification models with high- 
order interactions, with application to logistic sequence prediction models and to logistic 
classification models. We conclude separately at the end of each chapter. 




(1.11) 



Chapter 2 

Avoiding Bias from Feature Selection 



Abstract. For many classification and regression problems, a large number of features are 
available for possible use — this is typical of DNA microarray data on gene expression, for 
example. Often, for computational or other reasons, only a small subset of these features 
are selected for use in a model, based on some simple measure such as correlation with 
the response variable. This procedure may introduce an optimistic bias, however, in which 
the response variable appears to be more predictable than it actually is, because the high 
correlation of the selected features with the response may be partly or wholly due to chance. 
We show how this bias can be avoided when using a Bayesian model for the joint distribution 
of features and response. The crucial insight is that even if we forget the exact values of the 
unselected features, we should retain, and condition on, the knowledge that their correlation 
with the response was too small for them to be selected. In this paper we describe how this 
idea can be implemented for "naive Bayes" and mixture models of binary data. Experiments 
with simulated data confirm that this method avoids bias due to feature selection. We also 
apply the naive Bayes model to subsets of data relating gene expression to colon cancer, and 
find that correcting for bias from feature selection does improve predictive performance. 



Part of this Chapter appeared as a technical report coauthored with Jianguo Zhang and Radford Neal. 
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2.1 Introduction 

Regression and classification problems that have a large number of available "features" (also 
known as "inputs", "covariates" , or "predictor variables") are becoming increasingly com- 
mon. Such problems arise in many application areas. Data on the expression levels of tens 
of thousands of genes can now be obtained using DNA microarrays, and used for tasks such 
as classifying tumors. Document analysis may be based on counts of how often each word in 
a large dictionary occurs in each document. Commercial databases may contain hundreds 
of features describing each customer. 

Using all the features available is often infeasible. Using too many features can result in 
"overfitting" when simple statistical methods such as maximum likelihood are used, with the 
consequence that poor predictions are made for the response variable (e.g., the class) in new 
items. More sophisticated Bayesian methods can avoid such statistical problems, but using 
a large number of features may still be undesirable. We will focus primarily on situations 
where the computational cost of looking at all features is too burdensome. Another issue 
in some applications is that using a model that looks at all features will require measuring 
all these features when making predictions for future items, which may sometimes be costly. 
In some situations, models using few features may be preferred because they are easier to 
interpret. 

For the above reasons, modellers often use only a subset of features, chosen by some 
simple indicator of how useful they might be in predicting the response variable — see, for 
example, the papers in (Guyon, et al. 2006). For both regression problems with a real- valued 
response variable and classification problems with a binary (0/1) class variable, one suitable 
measure of how useful a feature may be is the sample correlation of the feature with the 
response. If the absolute value of this sample correlation is small, we might decide to omit 
the feature from our model. This criterion is not perfect, of course — it may result in a 
relevant feature being ignored if its relationship with the response is non-linear, and it may 
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result in many redundant features being retained even when they all contain essentially the 
same information. Sample correlation is easily computed, however, and hence is an attractive 
criterion for screening a large number of features. 

Unfortunately, a model that uses only a subset of features, selected based on their high 
correlation with the response, will be optimistically biased — i.e., predictions made using the 
model will (on average) be more confident than is actually warranted. For example, we might 
find that the model predicts that certain items belong to class 1 with probability 90%, when in 
fact only 70% of these items are in class 1. In a situation where the class is actually completely 
unpredictable from the features, a model using a subset of features that purely by chance had 
high sample correlation with the class may produce highly confident predictions that have 
less actual chance of being correct than just guessing the most common class. The feature 
selection bias has also been noticed in the literature by a few researchers, see for example, the 
papers (Ambroise and McLachlan 2002), (Lecocke and Hess 2004), (Singhi and Liu 2006), and 
(Raudys, Baumgartner and Somorjai 2005). They pointed out that if the feature selection 
is performed externally to the cross-validation assessment (ie, cross-validation is applied to 
a subset of features selected in advance based on all observations), the classification error 
rate will be highly underestimated (could be 0%). It is therefore suggested that feature 
selection should be performed internally to the cross-validation procedure, ie, re-selecting 
features whenever the training set and test set are changed. This modified cross-validation 
procedure avoids underestimating the error rate and assesses properly the predictive method 
plus the feature selection method. However, it does not provide a scheme for constructing 
a better predictive method that can give out well-calibrated predictive probabilities for test 
cases. We propose a Bayesian solution to this problem. 

This optimistic bias comes from ignoring a basic principle of Bayesian inference — that 
we should base our conclusions on probabilities that are conditional on all the available infor- 
mation. If we have an appropriate model, this principle would lead us to use all the features. 
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This would produce the best possible predictive performance. However, we assume here 
that computational or other pragmatic issues make using all features unattractive. When we 
therefore choose to "forget" some features, we can nevertheless still retain the information 
about how we selected the subset of features that we use in the model. Properly conditioning 
on this information when forming the posterior distribution eliminates the bias from feature 
selection, producing predictions that are as good as possible given the information in the 
selected features, without the overconfidence that comes from ignoring the feature selection 
process. 

We can use the information from feature selection procedure only when we model the 
features and the response jointly. We show in this Chapter this information can be easily 
incorporated into our inference in a Bayesian framework. We particularly apply this method 
to naive Bayes models and mixture models. 

2.2 Our Method for Avoiding Selection Bias 

Suppose we wish to predict a response variable, y, based on the information in the numerical 
features x±, . . . , x p , which we sometimes write as a vector, x. Our method is applicable both 
when y is a binary (0/1) class indicator, as is the case for the naive Bayes models discussed 
later, and when y is real-valued. We assume that we have complete data on n "training" 
cases, for which the responses are y^\ . . . , y( n > (collectively written as |/ train ) and the feature 
vectors are x^\ . . . , x^ (collectively written as x tI3 - in ). (Note that when y, x, or x t are 
used without a superscript, they will refer to some unspecified case.) We wish to predict 
the response for one or more "test" cases, for which we know only the feature vector. Our 
predictions will take the form of a distribution for y, rather than just a single-valued guess. 

We are interested in problems where the number of features, p, is quite big — perhaps 
as large as ten or a hundred thousand — and accordingly (for pragmatic reasons) we intend 
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to select a subset of features based on the absolute value of each feature's sample correlation 
with the response. The sample correlation of the response with feature t is defined as follows 
for as zero if the denominator below is zero): 



£ (y (l) - y) (x? ] - xt) 

ln l=1 7= (2.1) 

E {v®-v) a JxkP-St) 2 



where y — - E and x t = \ E x f' '• The numerator can be simplified to E (y^ —y)x[ l \ 

i=l i=l i=l 

Although our interest is only in predicting the response, we assume that we have a 
model for the joint distribution of the response together with all the features. From such 
a joint distribution, with probability or density function P(y, x%, . . . , x p ), we can obtain the 
conditional distribution for y given any subset of features, for instance P(y \ x\, . . . , x k ), with 
k < p. This is the distribution we need in order to make predictions based on this subset. 
Note that selecting a subset of features makes sense only when the omitted features can be 
regarded as random, with some well-defined distribution given the features that are retained, 
since such a distribution is essential for these predictions to be meaningful. This can be seen 
from the following expression: 

P(y | xi, . . . , x fc ) = J "J P (y \xi,...,x k , x k+ i, . . . , x p ) ■ 

P{x k+1 , . . . , x p I x x , . . . , x k ) dx k+ i ■■■dx p (2.2) 

If P(xfc + i, . . . , x p I Xi, . . . , x k ) does not exist in any meaningful sense — as would be the case, 
for example, if the data were collected by an experimenter who just decided arbitrarily what 
to set Xk+i, ■ ■ ■ , x p to — then P(y \ x±,..., x k ) will also have no meaning. 

Consequently, features that cannot usefully be regarded as random should always be 
retained. Our general method can accommodate such features, provided we use a model 
for the joint distribution of the response together with the random features, conditional on 



2 Avoiding Bias from Feature Selection 



18 



given values for the non-random features. However, for simplicity, we will ignore the possible 
presence of non-random features in this paper. 

We will assume that a subset of features is selected by fixing a threshold, 7, for the 
absolute value of the correlation of a selected feature with the response. We then omit 
feature t from the feature subset if |COR(y train , xj rain )| < 7, retaining those features with a 
greater degree of correlation. Another possible procedure is to fix the number of features, k, 
that we wish to retain, and then choose the k features whose correlation with the response 
is greatest in absolute value, breaking any tie at random. If s is the retained feature with 
the weakest correlation with the response, we can set 7 to |COR(y train , x* rain )|, and we will 
again know that if t is any omitted feature, |COR(y train , x£ rain )| < 7. If either the response 
or the features have continuous distributions, exact equality of sample correlations will have 
probability zero, and consequently this situation can be treated as equivalent to one in which 
we fixed 7 rather than k. If sample correlations for different features can be exactly equal, 
we should theoretically make use of the information that any possible tie was broken the 
way that it was, but ignoring this subtlety is unlikely to have any practical effect, since ties 
are still likely to be rare. 

Regardless of the exact procedure used to select features, we will denote the number of 
features retained by k, we will renumber the features so that the subset of retained features 
is 37, . . . , Xk, and we will assume we know that |COR(?/ train , Xj rain )| < 7 for t = k+1, . . . ,p. 

We can now state the basic principle behind our bias-avoidance method: When forming 
the posterior distribution for parameters of the model using a subset of features, we should 
condition not only on the values in the training set of the response and of the k features we 
retained, but also on the fact that the other p—k features have sample correlation with the 
response that is less than 7 in absolute value. That is, the posterior distribution should be 
conditional on the following information: 

y train , a^X, |COR(y train , x* rain )| < 7 for t = k+1, . . . ,p (2.3) 
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where x^.p = (x^ 111 , . . . , x^ rain ). 

We claim that this procedure of conditioning on the fact that selection occurred will 
eliminate the bias from feature selection. Here, "bias" does not refer to estimates for model 
parameters, but rather to our estimate of how well we can predict responses in test cases. 
Bias in this respect is referred to as a lack of "calibration" — that is, the predictive proba- 
bilities do not represent the actual chances of events (Dawid 1982). If the model describes 
the actual data generation mechanism, and the actual values of the model parameters are 
indeed randomly chosen according to our prior, Bayesian inference always produces well- 
calibrated results, on average with respect to the data and model parameters generated from 
the Bayesian model. The proof that the Bayesian inference is well-calibrated is given in the 
Appendix 1 to this Chapter. 

In justifying our claim that this procedure avoids selection bias (ie, is well-calibrated), 
we will assume that our model for the joint distribution of the response and all features, 
and the prior we chose for it, are appropriate for the problem, and that we would therefore 
not see bias if we predicted the response using all the features. Now, imagine that rather 
than selecting a subset of features ourselves, after seeing all the data, we instead set up 
an automatic mechanism to do so, providing it with the value of 7 to use as a threshold. 
This mechanism, which has access to all the data, will compute the sample correlations of 
all the features with the response, select the subset of features by comparing these sample 
correlations with 7, and then erase the values of the omitted features, delivering to us only the 
identities of the selected features and their values in the training cases. If we now condition 
on all the information that we know, but not on the information that was available to the 
selection mechanism but not to us, we will obtain unbiased inferences. The information we 
know is just that of (2.3) above. 

The class of models we will consider in detail may include a vector of latent variables, z, 
for each case. Model parameters 9±, . . . , 9 P (collectively denoted 6) are associated with the 
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t=l..p 




i = l..n 



Figure 2.1: A directed graphical model for the general class of models we are considering. 
Circles represent variables, parameters, or hyperparameters. Arrows represent possible direct 
dependencies (not all of which are necessarily present in all models in this class). The 
rectangles enclose objects that are repeated; an object in both rectangles is repeated in both 
dimensions. The case index, i, is shown as ranging over the n training cases, but test cases 
(not shown) belong in this rectangle as well. This diagram portrays a model where cases are 
independent given a and 6, though this is not essential. 

p features; other parameters or hyperparameters, a, not associated with particular features, 
may also be present. Conditional on and a, the different cases may be independent, 
though this is not essential for our method. Our method does rely on the values of different 
features (in all cases) being independent, conditional on 6, a, y train ; and z train . Also, in the 
prior distribution for the parameters, $i, . . . , 9 p are assumed to be conditionally independent 
given a. These conditional independence assumptions are depicted graphically in Figure 2.1. 

If we retain all features, our prediction for the response, y*, in a test case for which we 
know the features, x* = (xj, . . . , x*), can be found from the joint predictive distribution for 
y* and x* given the data for all training cases, written as y train and £c train : 



The posterior, P(a, 6 | y train ; jc train ), is proportional to the product of the prior and the like- 



P(y* I x*, y 



train 




P(y*, x* | y train , a; train ) 

P (^X* I ^ tra ^ n ^gtrain^ 

/ / P(y*, x* | a, 0) P(a, \ y train , as train ) dadO 



(2.4) 



(2.5) 



/ / P(x* | a, 6) P(a, 6 \ y train , cc train ) da dO 
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lihood: 

P(a, | ?/ train , a? train ) oc P(a, 0) P(y train , a; train | a, 6) (2.6) 

oc P{a) P(0* | a) J] P(y« | a, 6) (2.7) 

i=l i=l 

where the second expression makes use of the conditional independence properties of the 
model. 

When we use a subset of only k features, the predictive distribution for a test case will 

be 

P(y*\x* 1:k , xlT,S) 



Jfp(y*, 


e v.k) P(<*, Ol.k 


y train , xf.% 11 , S)dadO vk 






y 


tram ; jjjtram^ dad0 1:k 



(2.8) 



where 5 represents the information regarding selection from (2.3), namely | COR(y train , x^ rain ) | < 
7 for t = k+1, • • • ,p. The posterior distribution for a and 0\± needed for this prediction can 
be written as follows, in terms of an integral (or sum) over the values of the latent variables, 

,y train . 

(2.9) 
(2.10) 

a, z train , y train ) dz train (2.11) 
Here again, the conditional independence properties of the model justify removing the con- 
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ditioning on 9 1 . k and x 1 ™^ 1 in the last factor. 

Computation of P(S | a, 2 train , y train ) ; which adjusts the likelihood to account for feature 
selection, is crucial to applying our method. Two facts greatly ease this computation. First, 
the Xt are conditionally independent given a, z, and y, which allows us to write this as a 
product of factors pertaining to the various omitted features. Second, these factors are all 
the same, since nothing distinguishes one omitted feature from another. Accordingly, 

p 

P(<S | a, z tiain , y tia[n ) = J P(|COR(y train ,xr in )| < 7 | a, z train , y train ) (2.12) 

t=k+l 

P(|COR(y train ,xr in )l <l\a, 2 trah \ y train ) P k (2.13) 

where in the second expression, t represents any of the omitted features. Since the time 
needed to compute the adjustment factor does not depend on the number of omitted features, 
we may hope to save a large amount of computation time by omitting many features. 

Computing the single factor we do need is not trivial, however, since it involves integrals 
over 9t and xf ain . We can write 

P(|COR(y train ,x* rain )| < 7 | a, z train , y train ) 

= J P(6 t | a) P(|COR(y train , < ain )| < 7 1 a, 6 t , z train , y train ) d6 t (2.14) 

Devising ways of efficiently performing this integral over 9 t and the integral over x£ rain im- 
plicit in the probability statement occurring in the integrand will be the main topic of our 
discussion of specific models below. 

Once we have a way of computing this factor, we can use standard Markov chain Monte 
Carlo (MCMC) methods to sample from P(a, 1:k , 2 train | y train , a^", S). The resulting sam- 
ple of values for a and 6 vk can be used to make predictions using equation (2.8), by ap- 
proximating the integrals in the numerator and denominator by Monte Carlo estimates. For 
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the naive Bayes model we will discuss in Section 2.3, however, Monte Carlo methods are 
unnecessary — a combination of analytical integration and numerical quadrature is faster. 

2.3 Application to Bayesian Naive Bayes Models 

In this chapter we show how to apply the bias correction method to Bayesian naive Bayes 
models in which both the features and the response are binary. Binary features are natural 
for some problems (e.g., test answers that are either correct or incorrect), or may result 
from thresholding real-valued features. Such thresholding can sometimes be beneficial — in 
a document classification problem, for example, whether or not a word is used at all may 
be more relevant to the class of the document than how many times it is used. Naive Bayes 
models assume that features are independent given the response. This assumption is often 
incorrect, but such simple naive Bayes models have nevertheless been found to work well for 
many practical problems (see for example Li and Jain 1998, Vaithyanathan, Mao, and Dom 
2000, Eyheramendy, Lewis, and Madigan 2003). Here we show how to correct for selection 
bias in binary naive Bayes models, whose simplicity allows the required adjustment factor to 
be computed very quickly. Simulations reported in Section 2.3.5 show that substantial bias 
can be present with the uncorrected method, and that it is indeed corrected by conditioning 
on the fact that feature selection occurred. We then apply the method to real data on gene 
expression relating to colon cancer, and again find that our bias correction method improves 
predictions. 

2.3.1 Definition of the Binary Naive Bayes Models 

Let = (xy , • • • , Xp , ) be the vector of p binary features for case i, and let be the binary 
response for case i, indicating the class. For example, = 1 might indicates that cancer 
is present for patient i, and yW = Q indicate that cancer is not present. Cases are assumed 
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to be independent given the values of the model parameters (ie, exchangeable a priori). 
The probability that y = 1 in a case is given by the parameter ip. Conditional on the class 
y in some case (and on the model parameters), the features assumed to be 

independent, and to have Bernoulli distributions with parameters <j) y< i, ■ ■ ■ , y , p , collectively 
written as <p y , with <fi = (<p , 1 ) representing all such parameters. Figure 2.2 displays the 
models. Formally, the data is modeled as 

yW |^ ~ Bernoulli (ip), for i = 1, . . . , n (2-15) 
Xj | y^\ 4> ~ Bernoulli for i — 1, . . . , n and j — 1, . . . ,p (2.16) 

We use a hierarchical prior that expresses the possibility that some features may have 
almost the same distribution in the two classes. In detail, the prior has the following form: 
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i> ~ Beta (/i,/ ) (2.17) 

a ~ Inverse-Gamma(a, b) (2-18) 

01,..., P iid Uniform(0, 1) (2.19) 

0oj) I a ) 0j IID Beta(a0j, — 9j)), for j = l,...,p (2.20) 



The hyperparameters = (9±, . . . ,9 P ) are used to introduce dependence between O • 
and 0i j, with a controlling the degree of dependence. Features for which o ,j an d 4>ij differ 
greatly are more relevant to predicting the response. When a is small, the variance of the 
Beta distribution in (2.20), which is 9j (l—9j) / (a+1), is large, and many features are likely to 
have predictive power, whereas when a is large, it is likely that most features will be of little 
use in predicting the response, since 0o,j an d 4>i,j are likely to be almost equal. We chose 
an Inverse-Gamma prior for a (with density function proportional to a~^ 1+a ^ exp(— b/a)) 
because it has a heavy upward tail, allowing for the possibility that a is large. Our method 
of correcting selection bias will have the effect of modifying the likelihood in a way that 
favors larger values for a than would result from ignoring the effect of selection. 

2.3.2 Integrating Away ip and </> 

Although the above model is defined with ip and tf> parameters for better conceptual under- 
standing, computations are simplified by integrating them away analytically. 

Integrating away tp, the joint probability of y train = (y^\ . . . ,y( n >) is as follows, where 
J( • ) is the indicator function, equal to 1 if the enclosed condition is true and if it is false: 

1 n n 

P{y—) = / v Y f l f \ ^-W°^ h (i-#- # (2.21) 

Jo 1 UoJJ- \Ji) 

= u(f u /o, E i(y {l) = i), E i{y {l) = o)) (2.22) 
v i=i i=i 7 
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The function U is defined as 

U(f f n n) - r (/o + A) iU + noTO + m) 

£/(/!, /o,n 1; no) - fWTOruTTiT^x) (2 ' 23) 

no ni 

nc/o+i-i) nc/i+€-i) 

= ^ ^ (2.24) 

n (/o+A+«-i) 

The products above have the value one when the upper limits of n or n\ are zero. The joint 
probability of y train and the response, y*, for a test case is similar: 



, n n \ 

= u(f u fo, £ i{y {t) = i) + Hy* = i), £ i{v (i) = °) + = °) ) ( 2 - 25 ) 

v i=l i=l y 

Dividing p(y train ; y*) by p(y train ) gives 

P(y* | y train ) = Bernoulli (2.26) 
Here, Bernoulli (y; ^) = ^ v (1 - and $ = (/i + JVi) / (/„ + fi + n), with JV„ = 

n 

£ J(yW = y). Note that xj) is just the posterior mean of tp based on y^\ . . . , y( n \ 
t=i 

Similarly, integrating over 0o,j and </>ij, we find that 



P(xf in | 0,-, a, y train ) = J] U{ad v a(l-0,), O w -) (2.27) 

y=0 



where = £ I(y {i) = y, xf = 0) and I yJ = £ 7(y w = y, = 1). 

i=l i=l 

With ip and ^ integrated out, we need deal only with the remaining parameters, a and 
6. Note that after eliminating ip and the (f>, the cases are no longer independent (though 
they are exchangeable). However, conditional on the responses, y train ; and on a, the values of 
different features are still independent. This is crucial to the efficiency of the computations 
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described below. 



2.3.3 Predictions for Test Cases using Numerical Quadrature 

We first describe how to predict the class for a test case when we are either using all features, 
or using a subset of features without any attempt to correct for selection bias. We then 
consider how to make predictions using our method of correcting for selection bias. 

Suppose we wish to predict the response, y*, in a test case for which we know the retained 
features x* vk = (xl, ■ ■ ■ ,x* k ) (having renumbered features as necessary). For this, we need 
the following predictive probability: 

ply* I ytrain\ p ( X * I y* ^rain ytrain\ 

P(V* I x* 1:k , xlf 1 , y train ) = — LAJ±1 1± 1 (2.28) 

E p (y* = y I y train ) p ( x l k \y* = y> x Tk^ y train ) 

Ie, we evaluate the numerator above for y* = and y* = 1, then divide by the sum to 
obtain the predictive probabilities. The first factor in the numerator, P(y* |y train ), is given 
by equation (2.26). It is sufficient to obtain the second factor up to a proportionality constant 
that doesn't depend on y*, as follows: 

P(x* /£ train I y* ytrairA 

P( x t.k I y*, x it, 2/ train ) = Z ™ i ! , ( 2 - 29 ) 

L> ( ™tram n.train \ 

r \ x i: k i y I 
oc p{x\^ x\t I y\ y train ) (2-30) 

This can be computed by integrating over a, noting that conditional on a the features are 
independent: 

P(x*l:k> <:t I V*, V***) = / P ( X *l:ki X Xt I «, V* , V***) ^ ( 2 -31) 

= / p(<*) n p ( x h x T n I «. y*> y tiain ) da ( 2 - 32 ) 

J 3=1 
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Each factor in the product above is found by using equation (2.27) and integrating over 6y. 



P(x*, x} rain | a, y*, y train ) 

= / P{x* | 9j, a, xf*, y train , y*) P{xf in \ 6,, a, y train ) dQ 3 (2.33) 
Jo 

r 1 1 

= / Bernoulli (a:*; r j) ] [ U(a6 j7 a(l-0j), I yj , O yj )d9j (2.34) 
Jo y=o 

where (f) y *j = (a9j + I y *j) / (a + N y *), the posterior mean of <p y *j given a and 9j. 

When using k features selected from a larger number, p, the predictions above, which 
are conditional on only x^ 3 ^ and y train ; are not correct - - we should also condition on 
the event, <S, that |COR(y train , x^ rain )| < 7 for j = k + l,...,p. We need to modify 
the predictive probability of equation (2.28) by replacing P(x\. k \ y*, xf.^, y train ) with 
P(x* 1:k I y*, x^ k h \ y train , S), which is proportional to P(x* 1:k , xffi, S \ y*, ?/ train ). Analo- 
gously to equations (2.31) and (2.32), we obtain 

P(x* ltk , x%£, S I y\ y train ) 
= J P(a) P(x* 1:k , x%?, S I a, y\ y^) da (2.35) 

= / P(a) P(S I a, y^) J] P(x*, xf in \ a, y\ y^) da (2.36) 
J 3=1 

The factors for the k retained features are computed as before, using equation (2.34). The 
additional correction factor that is needed (presented earlier as equation (2.13)) is 

p 

P(S I a, y train ) = ' P(|COR(?/ train ,x} rain )| < 7 I a, y train ) (2.37) 
j=k+i 

= [p(|COR(y train ,xr in )| < 7 I a, y train ) P (2.38) 



where t is any of the omitted features, all of which have the same probability of having a 
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small correlation with y. We discuss how to compute this adjustment factor in the next 
section. 

To see intuitively why this adjustment factor will correct for selection bias, recall that as 
discussed in Section (2.3.1), when a is small, features will be more likely to have a strong 
relationship with the response. If the likelihood of a is based only on the selected features, 
which have shown high correlations with the response in the training dataset, it will favor 
values of a that are inappropriately small. Multiplying by the adjustment factor, which 
favors larger values for a, undoes this bias. 

We compute the integrals over a in equations (2.32) and (2.36) by numerical quadrature. 
We use the midpoint rule, applied to u — F(a), where F is the cumulative distribution 
function for the Inverse-Gamma (a, b) prior for a. The prior for u is uniform over (0, 1), and 
so needn't be explicitly included in the integrand. With K points for the midpoint rule, the 
effect is that we average the value of the integrand, without the prior factor, for values of a 
that are the 0.5/ K, 1.5/ K, ... ,1 — 0.5/ K quantiles of its Inverse-Gamma prior. For each a, 
we use Simpson's Rule to compute the one-dimensional integrals over 9j in equation (2.34). 



2.3.4 Computation of the Adjustment Factor for Naive Bayes 
Models 

Our remaining task is to compute the adjustment factor of equation (2.38), which de- 
pends on the probability that a feature will have correlation less than 7 in absolute value. 
Computing this seems difficult — we need to sum the probabilities of cc£ rain given y train ; a 
and $t over all configurations of xf ain for which |COR(y train , xj rain )| < 7 — but the com- 
putation can be simplified by noticing that COR(a^ rain , y train ) can be written in terms of 
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h = EI=i% W = 0, xf ] = 1) and h = £" =1 %« = h 4 ] = 1), as follows: 

n 

E (•'/" - y) ^ 

COR(xr in , y^) = 1=1 (2.39) 

(2.40) 



i=i v »=i 

(0 - 17) J + 



Vn^l-y) v / /o + /i-(/o + /i) 2 /^ 

We write the above as Cor(/ , h, y), taking n as known. This function is defined for < 
Iq < n(l— y) and < I± < ny. 



Fixing n, y, and 7, we can define the following sets of values for Iq and I\ (for some 



feature x t ) in terms of the resulting correlation with y: 

L = {(J ,/i) : Cor(J ,/ 1 ,27)=0} (2.41) 

L+ = {(I ,h) : 0<Cor(J ,Ii,y)<7} (2-42) 

L_ = {(10,70 = -7<Cor(/ ,/ 1 ,y)<0} (2.43) 

H + = {(/ ,/0 : 7<Cor(/ ,/i,17)} (2.44) 

= {(Jo,/!) : Cor(/ ,J 1 ,y)<- 7 } (2.45) 



A feature will be discarded if (J , ii) G L_ U L U L + and retained if (J , Ji) G if_ U //+. 
These sets are illustrated in Figure 2.3. 



We can write the probability needed in equation (2.38) using either L_, L , and L + or 
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Figure 2.3: The Cor function for a dataset with n = 22 and y = 14/22. The values of 
Cor(/ , h,y) are shown for the valid range of I and I\. Using 7 = 0.2, the values of (IoJi) 
in Lq are shown in dark grey, those in L_ or L + in medium grey, and those in i7_ or H + in 
light grey 

H_ and H + . We will take the latter approach here, as follows: 

P{ |COR(xr in ,Z/ train )l < 7 I «, y train ) 
= 1 - P((I ,h) e H„UH+ I a, y^) (2.46) 

= 1 - p ( J °> J i 1 a > ^ train ) ( 2 - 47 ) 

(io,/i)e-ff-uff + 

We can now exploit symmetries of the prior and of the Cor function to speed up compu- 
tation. First, note that Cor(/ , h,V) = — Cor(n(l— y) — 1 , ny — Ii,y), as can be derived from 
equation (2.40), or by simply noting that exchanging labels for the classes should change only 
the sign of the correlation. The one-to-one mapping (7 , h) — » (^(1— y) — lo, ny — h), which 
maps H_ and H + and vice versa (similarly for L_ and L + ), therefore leaves Cor unchanged. 
The priors for 9 and (see (2.19) and (2.20)) are symmetrical with respect to the class labels 
and 1, so the prior probability of (I , Ii) is the same as that of (n(l— y) — I , ny — Ii). We 
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can therefore rewrite equation (2.47) as 

P( |COR(^ rain ,y train )l < 7 I «, 2/ train ) = 1 ~ 2 Yl p ( J o ? J i I a > V^) ( 2 - 48 ) 

(/ ,Ji)efr+ 

At this point we write the probabilities for Iq and I\ in terms of an integral over t , and 
then swap the order of summation and integration, obtaining 

(I ,h)eH + J ° (I ,h)€H + 

The integral over 0$ can be approximated using some one-dimensional numerical quadrature 
method (we use Simpson's Rule), provided we can evaluate the integrand. 

The sum over H + can easily be delineated because Cor (Jo, Jl, V) is a monotonically de- 
creasing function of Iq, and a monotonically increasing function of Ji, as may be confirmed 
by differentiating with respect to Jo and I\. Let bo be the smallest value of I± for which 
Cot(0, I±,y) > 7. Taking the ceiling of the solution of Cor(0,Ji,I/) = 7, we find that 
bo = [l/(V n + (1 — V)/( n yi 2 ))] ■ F° r bo < h < ny, let be the largest value of Io for 
which Cor(J , h, y) > 7- We can write 

ny f/j 

£ P(J , ii I a, 9 t , y^) = J2 P ( /o ' h I a ' **> ^ tiain ) ( 2 - 5 °) 

(/o,Ii)Eff+ /i=6o /o=0 

Given a and ^, Jo an d A are independent, so we can reduce the computation needed by 
rewriting the above expression as follows: 
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Note that the inner sum can be updated from one value of I\ to the next by just adding any 
additional terms needed. This calculation therefore requires 1+ny — bo < n evaluations of 
P(Ii | a, 9 t , y train ) and l+r nV < n evaluations of P(I \ a, 9 t , y train ). 

To compute P(I\ \ a, 9t, y train ), we multiply the probability of any particular value for 
^.tram - m ^j,-.^ there are I\ cases with y = 1 and x t = 1 by the number of ways this can occur. 
The probabilities are found by integrating over o ,i an d 0i,t, as described in Section 2.3.2. 
The result is 

P(h | a, 6 t , y train ) = ( n y )u{a9 t , a(l-9 t ), h, ny - h) (2.52) 

Similarly, 

P(I | a, 9 U y^) = ( < l ~y) )u(a9 t , a(l-9 t ), J , n(l-y) - I ) (2.53) 

One can easily derive simple expressions for P(I\ \ a, 9 t , y train ) and P(Iq \ a, 9 t , y train ) in 
terms of P(h — 1 | a, 9 t , y train ) and P(Iq — 1 | a, 9 t , y train ) ; which avoid the need to compute 
gamma functions or large products for each value of Iq or I\ when these values are used 
sequentially, as in equation (2.51). 

2.3.5 A Simulation Experiment 

In this section, we use a dataset generated from the naive Bayes model defined in Section 2.3.1 
to demonstrate the lack of calibration that results when only a subset of features is used, 
without correcting for selection bias. We show that our bias-correction method eliminates 
this lack of calibration. We will also see that for the naive Bayes model only a small amount 
of extra computational time is needed to compute the adjustment factor needed by our 
method. 

Fixing a = 300, and p = 10000, we used equations (2.16), (2.19) and (2.20) to generate a 
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Figure 2.4: The absolute value of the sample correlation of each feature with the binary 
response, in the training set, and in the test set. Each dot represents one of the 10000 binary 
features. The training set correlations of the 1st, 10th, 100th, and 1000th most correlated 
features are marked by vertical lines. 
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set of 200 training cases and a set of 2000 test cases, both having equal numbers of cases with 
y = and y — 1. We then selected four subsets of features, containing 1, 10, 100, and 1000 
features, based on the absolute values of the sample correlations of the features with y. The 
smallest correlation (in absolute value) of a selected feature with the class was 0.36, 0.27, 
0.21, and 0.13 for these four subsets. These are the values of 7 used by the bias correction 
method when computing the adjustment factor of equation (2.38). Figure 2.4 shows the 
absolute value of the sample correlation in the training set of all 10000 features, plotted 
against the sample correlation in the test set. As can be seen, the high sample correlation 
of many selected features in the training set is partly or wholly a matter of chance, with the 
sample correlation in the test set (which is close to the real correlation) often being much 
less. The role of chance is further illustrated by the fact that the feature with highest sample 
correlation in the test set is not even in the top 1000 by sample correlation in the training 
set. 

For each number of selected features, we fit this data using the naive Bayes model with 
the prior for ip (equation (2.17)) having fo — fi — 1 and the prior for a (equation (2.18)) 
having shape parameter a = 0.5 and rate parameter b = 5. We then made predictions for the 
test cases using the methods described in Section 2.3.3. The "uncorrected" method, based on 
equation (2.28), makes no attempt to correct for the selection bias, whereas the "corrected" 
method, with the modification of equation (2.36), produces predictions that account for the 
procedure used to select the subset of features. We also made predictions using all 10000 
features, for which bias correction is unnecessary. 

We compared the predictive performance of the corrected method with the uncorrected 
method in several ways. First, we looked at the error rate when classifying test cases by 
thresholding the predictive probabilities at 1/2. As can be seen in Figure 2.5, there is little 
difference in the error rates with and without correction for bias. However, the methods 
differ drastically in terms of the expected error rate — the error rate we would expect based 
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number of features selected number of features selected 

Figure 2.5: Actual and expected error rates with varying numbers of features selected, with 
and without correction for selection bias. The solid line is the actual error rate on test cases. 
The dotted line is the error rate that would be expected based on the predictive probabilities. 
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Figure 2.6: Performance in terms of average minus log probability and average squared error, 
with varying numbers of features selected, with and without correction for selection bias. 
The left plot shows minus the average log probability of the correct class for test cases, with 
1, 10, 100, 1000, and all 10000 features selected. The dashed line is with bias correction, the 
dotted line without. The right plot is similar, but shows average squared error on test cases. 
Note that when all 10000 features are used, there is no difference between the corrected and 
uncorrected methods. 
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on the predictive probabilities for the test cases, equal to {l/N) £\ pW J(pW < 0.5) + 
(l-pW)J(pW > 0.5), where pW is the predictive probability of class 1 for test case z. The 
predictive probabilities produced by the uncorrected method would lead us to believe that we 
would have a much lower error rate than the actual performance. In contrast, the expected 
error rates based on the predictive probabilities produced using bias correction closely match 
the actual error rates. 

Two additional measures of predictive performance are shown in Figure 2.6. One measure 
of performance is minus the average log probability of the correct class in the N test cases, 
which is — (l/N) YliLi [v^ l°g(P^) + (1~ U^ 1 ) l°g(l — P^)]- This measure heavily penalizes 
test cases where the actual class has a predictive probability near zero. Another measure, 
less sensitive to such drastic errors, is the average squared error between the actual class (0 
or 1) and the probability of class 1, given by (l/N) Yli=i(v^ —p^) 2 - The corrected method 
outperforms the uncorrected one by both these measures, with the difference being greater for 
minus average log probability. Interestingly, performance of the uncorrected method actually 
gets worse when going from 1 feature to 10 features. This may be because the single feature 
with highest sample correlation with the response does have a strong relationship with the 
response (as may be likely in general), whereas some other of the top 10 features by sample 
correlation have little or no real relationship. 

We also looked in more detail at how well calibrated the predictive probabilities were. 
Table 2.1 shows the average predictive probability for class 1 and the actual fraction of cases 
in class 1 for test cases grouped according to the first decimal of their predictive probabilities, 
for both the uncorrected and the corrected method. Results are shown using subsets of 1, 
10, 100, and 1000 features, and using all features. We see that the uncorrected method 
produces overconfident predictive probabilities, either too close to zero or too close to one. 
The corrected method avoids such bias (the values for "Pred" and "Actual" are much closer), 
showing that it is well calibrated. 
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Table 2.1: Comparison of calibration for predictions found 
with and without correction for selection bias, on data simu- 
lated from the binary naive Bayes model. Results are shown 
with four subsets of features and with the complete data 
(for which no correction is necessary). The test cases were 
divided into 10 categories by the first decimal of the pre- 
dictive probability of class 1, which is indicated by the 1st 
column "C". The table shows the number of test cases in 
each category for each method ("#"), the average predic- 
tive probability of class 1 for cases in that category ( "Pred" ) , 
and the actual fraction of these cases that were in class 1 
("Actual"). 
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Figure 2.7: Posterior distributions of log(a) for the simulated data, with different numbers 
of features selected. The true value of log (a) is 5.7, shown by the vertical line. The solid line 
is the posterior density using all features. For each number of selected features, the dashed 
line is the posterior density including the factor that corrects for selection bias; the dotted 
line is the posterior density without bias correction. The dashed and solid lines overlap in 
the bottom two graphs. The dots mark the values of log (a) used to approximate the density, 
at the 0.5/K, 1.5/K, . . . , (K — 0.5) /K quantiles of the prior distribution (where K = 30). 
The probabilities of x train at each of these values for a were computed, rescaled to sum to K, 
and finally multiplied by the Jacobian, aP(a), to obtain the approximation to the posterior 
density of log (a) 



Number of Features Selected 1 10 100 1000 Complete data 



Uncorrected Method 11 19 107 1057 10639 

Corrected Method 12 19 107 1057 10639 



Table 2.2: Computation times from simulation experiments with naive Bayes models 
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The biased predictions of the uncorrected method result from an incorrect posterior 
distribution for a, as illustrated in Figure 2.7. Without bias correction, the posterior based 
on only the selected features incorrectly favours values of a smaller than the true value of 
300. Multiplying by the adjustment factor corrects this bias in the posterior distribution. 

Our software (available from http : //www.utstat .utoronto . ca/~longhai) is written in 
the R language, with some functions for intensive computations such as numerical integration 
and computation of the adjustment factor written in C for speed. We approximated the 
integral with respect to a using the midpoint rule with K = 30 values for F(a), as discussed 
at the end of Section 2.3.3. The integrals with respect to 9 in equations (2.34) and (2.49) 
were approximated using Simpson's Rule, evaluating 9 at 21 points. 

Computation times for each method (on a 1.2 GHz UltraSPARC III processor) are shown 
in Table 2.2. The corrected method is almost as fast as the uncorrected method, since the 
time to compute the adjustment factor is negligible compared to the time spent computing 
the integrals over 9j for the selected features. Accordingly, considerable time can be saved 
by selecting a subset of features, rather than using all of them, without introducing an 
optimistic bias, though some accuracy in predictions may of course be lost when we discard 
the information contained in the unselected features. 

2.3.6 A Test Using Gene Expression Data 

We also tested our method using a publicly available dataset on gene expression in normal 
and cancerous human colon tissue. This dataset contains the expression levels of 6500 genes 
in 40 cancerous and 22 normal colon tissues, measured using the Affymetrix technology. The 
dataset is available at http://geneexpression.cinj.org/~notterman/affyindex.html. 
We used only the 2000 genes with highest minimal intensity, as selected by Alon, Barkai, 
Notterman, Gish, Mack, and Levine (1999). In order to apply the binary naive Bayes model 
to the data, we transformed the real-value data into binary data by thresholding at the 
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median, separately for each feature. 

We divided these 2000 genes randomly into 10 equal groups, producing 10 smaller 
datasets, each with 200 features. We applies the corrected and uncorrected methods sepa- 
rately to each of these 10 datasets, allowing some assessment of variability when comparing 
performance. For each of these 10 datasets, we used leave-one-out cross validation to obtain 
the predictive probabilities over the 62 cases. In this cross-validation procedure, we left out 
each of the 62 cases in turn, selected the five features with the largest sample correlation 
with the response (in absolute value), and found the predictive probability for the left-out 
case using the binary naive Bayes model, with and without bias correction. The absolute 
value of the correlation of the last selected feature was around 0.5 in all cases. We used the 
same prior distribution, and the same computational methods, as for the demonstration in 
Section 2.3.5. 

Figure 2.8 plots the predictive probabilities of class 1 for all cases, with each of the 
10 subsets of features. The tendency of the uncorrected method to produce more extreme 
probabilities (closer to and 1) is clear. However, when the predictive probability is close to 
0.5, there is little difference between the corrected and uncorrected methods. Accordingly, the 
two methods almost always classify cases the same way, if prediction is made by thresholding 
the predictive probability at 0.5, and have very similar error rates. Note, however, that 
correcting for bias would have a substantial effect if cases were classified by thresholding 
the predictive probability at some value other than 0.5, as would be appropriate if the 
consequences of an error are different for the two classes. 

Figure 2.10 compares the two methods in terms of average minus log probability of the 
correct class and in terms of average squared error. From these plots it is clear that bias 
correction improves the predictive probabilities. In terms of average minus log probability, 
the corrected method is better for all 10 datasets, and in terms of average squared error, the 
corrected method is better for 8 out of 10 datasets. (A paired t test with these two measures 
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Figure 2.8: Scatterplots of the predictive probabilities of 
class 1 for the 10 subsets drawn from the colon cancer gene 
expression data, with and without correction for selection 
bias. Black circles are cases that are actually in class 1 
(cancer); hollow circles are cases that are actually in class 
0. Note that many case with predictive probabilities close 
to or 1 may overlap. 
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Figure 2.9: Actual versus expected error rates on the colon cancer datasets, with and without 
bias correction. Points are shown for each of the 10 subsets of features used for testing. 
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Figure 2.10: Scatterplots of the average minus log probability of the correct class and of the 
average squared error (assessed by cross validation) when using the 10 subsets of features 
for the colon cancer gene expression data, with and without correcting for selection bias. 
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produced p-values of 0.00007 and 0.019 respectively.) 

Finally, Figure 2.9 shows that our bias correction method reduces optimistic bias in 
the predictions. For each of the 10 datasets, this plot shows the actual error rate (in the 
leave-one-out cross-validation assessment) and the error rate expected from the predictive 
probabilities. For all ten datasets, the expected error rate with the uncorrected method is 
substantially less than the actual error rate. This optimistic bias is reduced in the corrected 
method, though it is not eliminated entirely. The remaining bias presumably results from 
the failure in this dataset of the naive Bayes assumption that features are independent within 

db clclSS. 

2.4 Application to Bayesian Mixture Models 

Mixture modelling is another way to model the joint distributions of the response variable and 
the predictor variables, from which we can find the conditional distribution of the response 
variable given the predictor variables. In this section we describe the application of the 
selection bias correction method to a class of binary mixture models, which is a generalization 
of the naive Bayes models. 

2.4.1 Definition of the Binary Mixture Models 

A complex distribution can be modeled using a mixture of finitely or infinitely many simple 
distributions, often called mixture components, for example, independent Gaussian distribu- 
tions for real values or independent Bernoulli distributions for binary values. Mixture models 
are often applied in density estimation, classification, and latent class analysis problems, as 
discussed, for example, by Everitt and Hand (1981), McLachlan and Basford (1988), and 
Titterington, et al. (1985). For finite mixture models with K components, the density or 
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z x o x i x P 






1 
1 
1 
1 

probability function of the observation x is written as 

K-l 

f(x \<f> ,..., <t> K -i) = Pk fk(x | 4> k ) (2.54) 

fc=0 

where p k is the mixing proportion of component and 4> k is the parameter associated with 
component fu- 

A Bayesian mixture model is often defined by introducing a latent label variable for each 
case, written as z. Given z = k, the conditional distribution of the observation x is the 
distribution for component k: 

x\ z=k,4> k ~ f k (x I <f) k ) (2.55) 

We consider a two-component mixture model in this section, which is a generalization 
of the naive Bayes model in Section 2.3. Most of the notation is therefore the same as in 
that section, except that we use the Oth feature, x , to represent the response y here (and so 
XQ ain is equivalent to y train ) for convenience of presentation. The parameters of the Bernoulli 
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distributions are (p Z) o (for the response Xq), and <f> Zt i, ■ ■ ■ ,<fi z ,p (for the features), collectively 
denoted as (f> z , for z = 0,1. The priors for O and 4> l are assigned in the same way as for 
binary naive Bayes model. The labels of the cases are assigned a prior in the same way as 
for y in naive Bayes models. Conditional on the component labels, all the features and the 
response are assumed to be mutually independent. The models are displayed by Figure 2.11, 
and are described formally as follows: 











Beta(/i,/ ) 


(2.56) 






z (n) | V 


IID 


Bernoulli (^) 


(2.57) 










Inverse-Gamma (ao, &o) 


(2.58) 






a 




Inverse-Gamma (a, 6) 


(2.59) 










r «o if J = 
a if j > 


(2.60) 






) ' " " ; 0p 


IID 


Uniform(0, 1) 


(2.61) 


<f>0,j,(f 






IID 


Beta (aij9j, aj(l — 6j)) j = 0, 1, • • • ,p 


(2.62) 


W i 








Bernoulli ((f) z (t)j) i = 1, • • • , n 


(2.63) 



The naive Bayes models described in Section 2.3 can be seen as a simpler form of the 
above mixture models by letting z^\ • ■ ■ , zW equal the responses :r rain , i.e., fixing 0o,o — 1 
and = 0. As for Bayesian naive Bayes models, ip and 9q, . . . ,9 p can be integrated away 
analytically from Bayesian binary mixture models. The cases are then no longer independent, 
but still are exchangeable. This integration reduces the number of the parameters, therefore 
improves Markov chain sampling or numerical quadrature if they are needed, though the 
resulting model may be harder to manipulate. 
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2.4.2 Predictions for Test Cases using MCMC 

Let us start with no attempt to correct for the selection bias. We want to predict the 
response, Xq, of a test case for which we know the retained features x\, ■ ■ ■ , x£ (renumbering 
the features as necessary), based on the training data x^ k in . For this, we need to calculate 
the predictive distribution: 

P{rr.* — 1 ™* I strain \ 

P(x* = 1 I x{. k , a#f ) = — l ^°~ 7 °- k) , , (2.64) 

V U I l.fc) U.fe J P(nf* _ 1 ~* ™train\ _|_ Df™* _ fl T>* strain \ V > 

JT{X — i, j, 1:k | x . fc ; -|- r i^x — u, «t, 1:fc | x Q . fc ) 

We therefore need to calculate P(x* Q . k \ x^ k ) for x* = 1 and Xq = 0. P{xl ± \ x^.p) can 
be written as: 



^train J OL J OL J 

P(0 0:kl ato, a, z train | a&f) d0 O: fc da da (2.65) 

^A^Orfc J strain ^ "0 ^ £» 

P«f | O:fcl a 0l q, Z trai ") P(e 0:fc ) P(«„) P(a) P( 2 trai ") 4, da da (2.66) 



The above integral is intractable analytically. We first approximate the integrals with 
respect to «o and a with the midpoint rule applied to the transformed variables uo = F (ao) 
and u = F(a), where F and F are the cumulative distribution functions of the priors for 
and a. Accordingly, the priors for uq and u are uniform over (0, 1). Suppose the midpoint 
rule evaluates the integrands with respect to Uq and u at K points respectively (K can be 
different for w an d u, for simplicity of presentation assume the same). After rewriting the 
summation with respect to uq and u over K points, (1 — 0.5)/ K, (2 — 0.5)/K, . . . , (K — 0.5)/K, 
in terms of cto and a, the midpoint rule approximation is equivalent to summing the integrand 
in (2.66), without the prior distribution for a and a, over the quantiles of the priors for 
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ao and a corresponding to probabilities (1 — 0.5) /K, (2 — 0.5) /K, . . . , (K — 0.5)/ K. Let us 
denote the K quantiles of ao by Ao, and denote the K quantiles of a by A. The integral 
in (2.66), with 1/P{xq±) omitted (since it is the same proportionality factor for all Xq), is 
approximated by: 

J2 E E / i <f,^,«o,«,2 tram )- 

P«f I Oo-.k, ao, a, z^) P(0 O:fc ) P(* train ) d0 O:fc (2.67) 

Since there are no prior terms in (2.67) for a an d a, the above approximation with 
midpoint rule applied to uq and u can also be seen as approximating the continuous Inverse- 
Gamma priors for ao and a by the uniform distributions over the finite sets Ao and A. 
Based on these discretized priors for a and a, we use Gibbs sampling method to draw 
samples from the posterior distribution of Oo-.k, ao, a, 2 train , allowing the integral in (2.67) to 
be approximated with the Monte Carlo method. The reason we use such a discretization for 
the prior for a is to ease the computation of the adjustment factor, which depends on a. As 
will be discussed later, when a is discrete, we can cache the values of the adjustment factors 
for future use when the same a is used again. 

We now start to derive the necessary formulae for performing Gibbs sampling for esti- 
mating (2.67). Using the results from Section 2.3.2, P{x* . k \ x^, 6o±, a, 2 train ) in (2.67) can 
be calculated as follows: 

P(x* . k | Oo-.k, a , a, z train ) 

i 

= P ( z * I z " ain ) P ( x o-k I <f , O :fc, «o, ot, Z^\ z*) (2.68) 

z*=0 

1 k 
= Y Bernoulli (z*; JJ Bernoulli (x*; 4> z * d ) (2.69) 

z*=0 j=0 

where ^ = (a^ + I^) / (a^+n^), j> = (/ 1 + n W)/(/ + /i + n), = £? =1 J(* = z), 
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I Z ,j = EtiH^ = z,xf = 1), and 0, d = EtiH^ = z,xf = 0). 

Again, using the results from Section 2.3.2, the distribution of given other training 
cases, denoted by which is needed to update with Gibbs sampling, can be found: 

k 

P(x^ k \x { -^,z^,e .. k ,ao,a) = JJ Bernoulli (xf;^.) (2.70) 

3=0 

where ftj. = {aft + + n^H) - 1), 7$. = £^ =1 = 1, z« = + i) 

and n^-V = £» =1 = s ± 

Similarly, the distribution of the whole training data given 2 train , 0:k , a o, a based on the 
k retained features can be found: 

k 1 

P«f | z*-, a ,a, 6 0:k ) = nil Ufa 6 * "iC 1 "^). 7 *>i> ( 2 - 71 ) 

j=0 2=0 

Integrating away ^ gives the prior for , • • ■ , z^ : 

P(z« • • • , * (n) ) = tf(/ 1) / 0) nM nM) (2.72) 

From the priors for zW, • ■ ■ , z^ n \ the conditional distribution of given all other z^ except 
zW, written as z^~ % \ can be found: 

P(z« | z (-0) = Bernoulli (z (i) ; (2.73) 

where ^H) = j^jgg and nW(-*) = £"=i 7(z s = 1, a ^ <) 

We now can write out the conditional distributions needed for performing Gibbs sampling. 
The conditional distribution of is proportional to the product of (2.73) and (2.70): 

k 

P(z® | a; train , z { -*\ e 0:k , ao, a) oc Bernoulli (z®; ] J Bernoulli {xf- $ { ~^.) (2.74) 

i=o 
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The conditional distribution of 6j is related to only feature j since the features are indepen- 
dent given z^ , ■ ■ • , z^ : 

1 

P{0j | sp.a^z'^) ex JJC/"(a0j-, a(l-8j), I zJ , O z J) (2.75) 

The conditional distribution of a given cc^ 11 , z tram , and #o:fc is proportional to the products 
of the factors for j > in (2.71) since the prior for a is uniform over A. And the conditional 
distribution of ao is proportional to the factor for j = in (2.71). 

The prediction described above is, however, invalid if the k features are selected from 
a large number. It needs to be modified to condition also on the information S, i.e., we 
should compute P{xq \ x*. k , Xq.^, S). The calculations are similar to the above, but with 
P(x%£> | O: fc, a, 2 train ) replaced by P(x^ n ,S | 0:kl a, z train ) in (2.66). Accordingly, the 
conditional distributions of a and z^\ ■ ■ ■ , z^ are multiplied by the following adjustment 
factor: 

P(S I x rain ,z train ,a) = (J P( |COR(xr in ,x rain )| < 7 | x rain , z train , 9 U a)dd t \ (2.76) 

Compared with the adjustment factor for the Bayesian naive Bayes models, the ad- 
justment factor (2.76) is more difficult to calculate, as we will discuss in the next section. 
Furthermore, this adjustment factor depends on both a and the unknown latent label vari- 
ables z^ l \ • • • , z( n \ for which we need to sample using Markov chain sampling method. We 
therefore need to recompute the adjustment factor whenever we change z^ l \ • • • , z^ during 
Markov chain sampling run. But we still need only to calculate the probability of one feature 
being discarded then raise it to the power of p — k. 



2 Avoiding Bias from Feature Selection 



51 



' 1 
1 
1 
1 



1 
1 

> 1 

Figure 2.12: Notations used in deriving the adjustment factor of Bayesian mixture models 
2.4.3 Computation of the Adjustment Factor for Mixture Models 

Computing P(S \ similar to what was done for Bayesian naive Bayes models 

in Section 2.3.4, with the difference that we condition on both £g ain and z train . The region 
|COR(xf ain , XQ ain )| < 7 can still be seen from Figure 2.3. The P(S | XQ ain , 2 train , a) is equal to 
the sum of P(Io, I\ | £Q ain , 2 train , a) over L + U L_ U L , or equivalently 1 minus the sum over 
H + U H^. The probability over H + is equal to the probability over if_ since the prior for 
9 t is symmetrical about 1/2. We therefore need to compute the probability for each point 
only in either H + or H_. We then exchange the summation over H + with the integration 
with respect to 9 t . Next we discuss only how to calculate P(Iq, h | Xg ain , z train , Q t , a) for each 
point (J , 7i). 

We divide the training cases according to 2 train into two groups, and let Iq = Y^h=i I(z® — 
z,Xq = 0,xf' = 1), and l{ 2 ' = I(z® = z,Xq = l,x[ 1 ' = 1), where z = 0,1. The probability 
of (Jq Z ',/{ z ') is found by summing over all configurations of feature t that have and 
results in 

P(I [ z \l[ z] | x^\z^\e u a) 

= ( S ){% )u(a9 t ,a(l-9 t ),lP + lP,nW-(lM + I?) (2.77) 
-'o *i 
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where riff 1 = ^2™ =l I(z^ = z,x^ = 0),n^ = Y17=i = z i x t = 1) an d n ^ = 

Then, the joint probability function of (Jo, h) is found by summing over all possible 
combinations of (if 1 , if) and (if 1 , if 1 ) that result in if + if 1 = i , if 1 + ij 1] = i x : 

l 

P(i , i! | < ain , 2 trah \ t , a) = ^ [J P (4 Z] > 4 Z] I < ah \ * trah \ 0*, ") (2.78) 

The way of finding the combinations of (if, if) and (if , if ) that satisfy if 1 +iQ 11 = io 
and if + if = ii is given in the Appendix 2 to this Chapter. 



2.4.4 A Simulation Experiment 

We tested our method using a data with 200 training cases and 2000 test cases, which are 
generated from a Bayesian mixture model, by setting a = 300, O o — 0.1, and 4>w = 0.9, and 
letting the number of z — 1 and z = be equal in both training and test sets. 

We then selected four subsets of features, containing 1, 10, 100, and 1000 features, based 
on the absolute values of the sample correlations of the features with y. The smallest 
correlation (in absolute value) of a selected feature with the class was 0.30, 0.24, 0.18, and 
0.12 for these four subsets. These are the values of 7 used by the bias correction method 
when computing the adjustment factor. 

For each number of selected features, we fit this data using the Bayesian mixture model 
with the prior for ip (equation (2.56)) having fo — fx — 1 and the Inverse-Gamma prior 
for both «o and a (equation (2.58) and (2.59)) both having shape parameter a = 0.5 and 
rate parameter 6 = 5. After using Gibbs sampling to train the model, with and without 
correction for the selection bias, we made predictions for the test cases. 

We compared the predictive performance of the methods with and without correction for 
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1 feature selected out of 10000 

Corrected Uncorrected 
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1083 


0.488 


0.472 


1083 


0.424 
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917 


0.536 


0.523 
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917 


0.617 
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0.472 
0.523 



10 features selected out of 10000 





Corrected 


Uncorrected 


# 


Pred 


Actual 


# 


Pred 


Actual 









19 


0.089 


0.105 









340 


0.155 


0.415 


19 


0.259 


0.105 


45 


0.269 


0.622 


317 


0.368 


0.420 


406 


0.359 


0.480 


530 


0.466 


0.494 


52 


0.424 


0.558 


552 


0.549 


0.500 


54 


0.560 


0.407 


480 


0.639 


0.529 


443 


0.649 


0.519 


100 


0.735 


0.620 


49 


0.735 


0.449 


2 


0.832 


1.000 


329 


0.854 


0.505 









263 


0.945 


0.593 



100 features selected out of 10000 



C 


1 

2 
3 
4 
5 
6 
7 



# 

71 
195 
237 
229 
234 
259 
253 
251 
192 

79 



Corrected 
Pred Actual 



0.072 
0.154 
0.250 
0.349 
0.454 
0.549 
0.652 
0.748 
0.848 
0.928 



0.183 
0.308 
0.300 
0.328 
0.504 
0.556 
0.565 
0.673 
0.729 
0.734 



# 

605 
122 
86 
63 
62 
90 
66 
87 
140 
679 



Uncorrected 
Pred Actual 



0.033 
0.144 
0.246 
0.350 
0.448 
0.551 
0.650 
0.749 
0.856 
0.965 



0.286 
0.361 
0.465 
0.508 
0.597 
0.589 
0.530 
0.529 
0.564 
0.666 



1000 features selected out of 10000 



# 
692 
129 
88 
64 
68 
71 
69 
S3 
143 
593 



Corrected 
Pred Actual 



Uncorrected 



0.033 
0.148 
0.247 
0.350 
0.454 
0.548 
0.646 
0.744 
0.857 
0.966 



0.140 
0.271 
0.375 
0.500 
0.426 
0.634 
0.580 
0.795 
0.804 
0.841 



# 

919 
33 
28 
21 
20 
25 
20 
31 
44 

859 



Pred Actual 



0.016 
0.145 
0.240 
0.350 
0.441 
0.552 
0.648 
0.749 
0.856 
0.980 



0.185 
0.455 
0.429 
0.619 
0.400 
0.480 
0.700 
0.645 
0.545 
0.818 



Table 2.3: Comparison of calibration for predictions found with and without correction for 
selection bias, on data simulated from a binary mixture model. The test cases were divided 
into 10 categories by the first decimal of the predictive probability of class 1, which is 
indicated by the 1st column "C" . The table shows the number of test cases in each category 
for each method ("#"), the average predictive probability of class 1 for cases in that category 
( "Pred" ) , and the actual fraction of these cases that were in class 1 ( "Actual" ) . 
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Figure 2.13: Actual and expected error rates with varying numbers (in log scale) of features 
selected, with and without correction for selection bias. The solid line is the actual error 
rate on test cases. The dotted line is the error rate that would be expected based on the 
predictive probabilities. 



Minus Average Log Probability Average Squared Error 




i 1 — i 1 — i 1 — r n 1 — i 1 — i 1 — r 
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Figure 2.14: Performance in terms of average minus log probability and average squared 
error, with varying numbers (in log scale) of features selected, with and without correction 
for selection bias. The left plot shows minus the average log probability of the correct class 
for test cases, with 1, 10, 100, and 1000 features selected. The dashed line is with bias 
correction, the dotted line without. The right plot is similar, but shows average squared 
error on test cases. 
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selection bias in several ways. Table 2.3 shows how well calibrated the predictive probabilities 
were, by looking at the actual fraction of class 1 of the test cases with predictive probabilities 
within each of the ten intervals evenly spaced in (0, 1). This shows that the methods with 
correction for selection bias are better calibrated than without correction. For the methods 
with correction for selection bias, the actual fractions are closer to predictive probabilities, 
whereas for the methods without such correction, they are more different, with predictive 
probabilities incorrectly close to or 1 for many cases. But we have seen that some bias, 
although less severe, still exists for the methods with correction, which we will explain later. 

The calibration can also be illustrated by comparing the actual error rate, from making 
predictions by thresholding the predictive probabilities at 0.5, to the expected error rate, 
equal to (1/iV) EiP (i)/ (P W < °- 5 ) + i 1 ~ p w )J(p (<) > 0.5), where is the predictive 
probability for test case i. As shown by Figure 2.13, the expected error rates for the methods 
without correction for selection bias are much lower than the actual error rates, showing that 
the predictions are overconfident. In contrast, the expected error rates and the actual error 
rates are much closer for the methods with correction for selection bias, though there are 
still gaps between them. From Figure 2.13, we also see that the actual error rates for the 
methods with and without correction for selection bias are almost the same. 

The remaining bias for the methods with correction for selection bias presumably results 
from Markov chain Monte Carlo method. Since the two groups are not very apart with 
a = 300, the latent values ■ ■ ■ , temporarily converge to the responses . . . , 
even with correction for selection bias. The estimates of </>oo and 0io are therefore very close 
to and 1. But the traces of a with correction for selection bias still move around much 
bigger values in A than without correction. The probabilities of a test case belonging to 
two groups with correction are therefore closer than without correction, making it difficult 
to decide the label of the test case. The correction method, implemented with simple Gibbs 
sampling, reduces the selection bias, but does not eliminate it entirely. This remaining bias 
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will be eliminated entirely if one uses a more sophisticated Markov chain sampling method 
that allows the Markov chains to explore more thoroughly in the space of z^\ ■ • ■ , z^ n \ Note 
that, however, this is a matter of computation rather than of theory. In theory, the bias 
will be eliminated entirely by conditioning on all information available in making Bayesian 
inference if the data sets are generated from the Bayesian model. 

The methods with and without correction for selection bias are also compared in terms 
of average minus log probability and average squared error, as shown in Figure 2.14. In 
both measures, the methods with the selection bias corrected are superior over the methods 
without correction. 

The predictive performance using the complete data was not shown in previous table 
and figures, because it is ironically worse than using selected features. This is because the 
Markov chain, using the simple Gibbs sampling, converges temporarily to only one group 
if the initial labels are drawn randomly from the permutations of 1, ... ,n. Again, a more 
sophisticated Markov chain sampling method will solve this problem. 

We used Gibbs sampling to train the model. In each iteration of Gibbs sampling, we 
update a and «o once, and repeat 5 time the combination of updating the labels z^ l \ ■ ■ ■ , 
once and updating each of Q\.^ 20 times. In approximating the continuous priors for a an d 
a with the uniform distribution over the quantiles of the priors (equation (2.67)), we chose 
K = 10, giving A = A = {2.60, 4.83, 7.56, 11.45, 17.52, 27.99, 48.57, 98.49, 279.60, 2543.14}. 
Simpson's Rule, which is used to approximate the integral with respect to 9 t for computing 
the adjustment factor, evaluates the integrand at 11 points. 

Our software (available from http://www.utstat.utoronto.ca/~longhai) is written 
entirely in R language. Computation times for each method (on a 2.2 GHz Opteron processor, 
running 50 iterations of Gibbs sampling as described above) are shown in Table 2.4. The 
computation of adjustment factor takes a large amount of extra time. This is because the 
computation of a single adjustment factor is more complex than naive Bayes models and 
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Number of Features Selected 



1 



10 



100 1000 



Uncorrected Method 



70 



82 



193 1277 



Corrected Method 



3324 2648 2881 4300 



Table 2.4: Computation times from simulation experiments with mixture models. 

the computation needs to be redone whenever the latent values z^ 1 ' , • • • , change. The 
current method for computing the adjustment factor can still be improved. However, the 
methods using selecting features and correcting for the selection bias still work faster than 
using complete data, which takes about 40000 seconds for updating 50 iterations. 

2.5 Conclusion and Discussion 

We have proposed a Bayesian method for making well-calibrated predictions for a response 
variable when using a subset of features selected from a larger number based on some measure 
of dependency between the feature and the response. Our method results from applying the 
basic principle that predictive probabilities should be conditional on all available information 
- in this case, including the information that some features were discarded because they 
appear weakly related to the response variable. This information can only be utilized when 
using a model for the joint distribution of the response and the features, even though we are 
interested only in the conditional distribution of the response given the features. 

We applied this method to naive Bayes models with binary features that are assumed to 
be independent conditional on the value of the binary response (class) variable. With these 
models, we can compute the adjustment factor needed to correct for selection bias. Crucially, 
we need only compute the probability that a single feature will exhibit low correlation with 
the response, and then raise this probability to the number of discarded features. Due to 
the simplicity of naive Bayes models, the methods with the selection bias corrected work 
as fast as the methods without considering this corrrection. Substantial computation time 
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can therefore be saved by discarding features that appear to have little relationship with the 
response. 

We also applied this method to mixture models for binary data. The computation of the 
adjustment factor is more complex than for naive Bayes models, and it needs to be computed 
many times. But the method is still feasible, and will be faster than using all features when 
the number of available features is huge. 

The practical utility of the bias correction method we describe would be much improved 
if methods for more efficiently computing the required adjustment factor could be found, 
which could be applied to a wide class of models. 

Appendix 1: 

Proof of the well- calibration of the Bayesian Prediction 

Suppose we are interested in predicting whether a random vector Y is in a set A if we know 
the value of another random vector X. Here, X is all the information we know for predicting 
Y, such as the information from the training data and the feature values of a test case. And 
Y could be any unknown quantity, for example a model parameters or the unknown response 
of a test case. For discrete Y, A may contain only a single value; for continuous Y, it is 
a set such that the probability of Y G A is not (otherwise any predictive method giving 
predictive probability is well-calibrated). From a Bayesian model for X and Y, we can 
derive a marginal joint distribution for X and Y, P(X, Y) (which may be a probability 
function or a density function, or a combination of probability and density function), by 
integrating over the prior for the model parameters. 

Let us denote a series of independent experiments from P(X,Y) as (Xi,Yi), for % = 
1,2, ... . Suppose a predictive method predicts that event Y e A will occur with probability 
Y(x) after seeing X = x. Y(x) is said to be well-calibrated if, for any two numbers 
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Ci, c 2 G (0, 1) (assuming c\ < c 2 ) such that P{ Y(Xi) G (c±, c 2 ) ) ^ 0, the fraction of Y i G A 
among those experiments with predictive probability, Y(Xi), between c\ and c 2 , will be 
equal to the average of the predictive proabilities (with P-probability 1), when the number 
of experiments, k, goes to oo, that is, 



This definition of well-calibration is a special case for iid experiments of what is defined 
in (Dawid 1982). Note that this concept of calibration is with respect to averaging over both 
the data and the parameters drawn from the prior. 

We will show that under the above definition of calibration, the Bayesian predictive 
function Y(x) = P(Y G A \ X = x) is well-calibrated. 

First, from the strong law of large numbers, the left-hand of (2.79) converges to: 



We then need only show that the expression (2.80) is actually equal to 0, i.e., the numerators 
in two terms are the same. This equality can be shown as follows: 



Eti IjYjEA and YpQ G (c 1; c 2 ) ) Eti fVQ H *W g fa, c 2 ) ) 
Eti /( Y(X t ) G ( Cl , c 2 ) ) Eti H Y(X t ) G ( Cl , c 2 ) ) 



(2.79) 



P(Y G A and Y{X) G (ci, c 2 )) E( F(X) J(F(X) G (ci, c 2 ) ) ) 
P(F(X) G ( Cl ,c 2 )) ^(^(X) G (d jC2 )) 



(2.80) 



P(Y eA&nd Y(X) G (ci,c 2 )) 




(2.81) 



(2.82) 



P(F(X)/(F(X)G( Cl ,c 2 ))) 



(2.83) 



What is essential from (2.81) to (2.82) is that the Bayesian predictive function Y(x) is just 



the conditional probability P(Y G A \ X = x). 



The Bayesian predictive function Y(x) = P(Y G A \ X = x) also has the following 
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property, which is helpful in understanding the concept of well-calibration: 

P(Y G A | Y(X) G (ci, c 2 )) = E(Y(X) \ Y(X) G (c 1; c 2 )) G (ci, c 2 ) (2.84) 

P(V G .4 | Y(X) G (ci, c 2 )) is just the first term in (2.80), and is equal to the second term 
in (2.80), which can be written as E(Y(X) \ Y{X) G (ci,c 2 )). This conditional expectation 
is obviously between c\ and c 2 . 

Appendix 2: 

Details of the Computation of the Adjustment Factor 
for Binary Mixture Models 

Delineating H + 

With the monotonicity of Cor(/ 1 , 7 ; XQ ain ) with respect to either I\ or I , we can easily 
determine the bound of I\ for each I that satisfies Cor(ii, I ] x) < 7 by solving the equation: 

|Cor(/ ,/ i; xr n )| =7 (2.85) 

then rounding to the appropriate integers and truncating them by and n\. I.e. the lower 
bound is max(0, |7] ) and the upper bound is min(ni, \u\), where I and u are the two solutions 
of equation (2.85). When J = 0, if I\ is also 0, we assume the correlation of £Q rain and xj rain 
is 0, therefore the lower bound is set to be 0. Similarly, when Jo = no, the upper bound is 
set to be ri\. 
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Computation of P(I^\I[ l | x rain , z train , 9 t , a) with formula (2.77) 

We need to calculate P(I [ Z \ l[ z] | x rain , z train , 6 h a) for all e {0, • ■ • ,n [ z] } and if 1 G 
{0, ■ • • , n^y. These values are saved with a matrix T"M for convenience. We don't have 
to evaluate each element of T"M by noting that the third factor of equation (2.78) depends 
only on 1^ + l[ z \ i.e. the number of x£ raln = 1. For each k G {0, 1, • • • , nffl + n^}, we 
need to evaluate it only once. Then go along the diagonal line Iq + I[ = k to obtain 
the elements of T^. The lower bound of /J 2 ' on this line is max(0, k — n^) and the upper 
bound is min(A;,n^). For each ij z ' between the lower and upper bounds, correspondingly 
if] =k- jW. 

For each line associated with k, only when jj 2 ' = max(0, k — nffl) we need to evaluate 
(2.78), then we can obtain the remaining values using the following relation between two 
successive elements on the line: 



tin \ l n 



P(li z \l[ z] | x rain ,z train ,a,^) V J o /V A 



P(iJ* J + 1, If - 1 | 4 rain , z train , a, t 



J Z| + 1 / V l[ z] - 1 
(^ + l)(nW-/W + l) 



(2.86) 



(2.87) 



Computation of P(J , Ji | x rain , z train , 6> 4 , a) with formula (2.78) 

For each (Jo, Ii) G -£/+, we need find all the pairs of {I^\ ij°') and (lj)\ 1^) that satisfy 

J^+Jo 11 = I ,and0<4 0l <n™ 0<4 1] <n 1] (2.88) 
= I uaa dO<I?<r^,0<I?<n? (2.89) 

The decompositions of Jo and I\ are independent and the methods are identical. Taking 
Iq as example, we determine the bound of Iq and 1^ by truncating the straight line of 
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7q + Iq — i with the square determined by (0, ) x (0, n [ ' ). By this way, we obtain the 
decompositions as follows: 

4°' e {max(0, J - ■ ■ • , min(4 0] , J )} = and if 1 = i - if 1 (2.90) 
i| 0] G {max(0, Ji - • • • , mm(nf ] , h)} = fig 1 , and = ii - if (2.91) 

We make a sub-matrix S*' ' from T^, which has been computed in Section 2.5, by taking 
the rows in B^} and the columns in B^, and accordingly make from T"M by taking 
the rows in Iq — B l S and the columns in ii — Bf}. Then multiplying and S 1 ' 1 ' element 
by element (i.e. the corresponding elements of and are multiplied together) makes 
matrix S. Summing all the elements of S together yields P(Iq, I\ | XQ ain , 2 train , 9 t , a). 



Chapter 3 

Compressing Parameters in Bayesian 
Models with High-order Interactions 

Abstract. Bayesian regression and classification with high order interactions is largely 
infeasible because Markov chain Monte Carlo (MCMC) would need to be applied with a huge 
number of parameters, which typically increases exponentially with the order. In this chapter 
we show how to make it feasible by effectively reducing the number of parameters, exploiting 
the fact that many interactions have the same values for all training cases. Our method uses 
a single "compressed" parameter to represent the sum of all parameters associated with 
a set of patterns that have the same value for all training cases. Using symmetric stable 
distributions as the priors of the original parameters, we can easily find the priors of these 
compressed parameters. We therefore need to deal only with a much smaller number of 
compressed parameters when training the model with MCMC. The number of compressed 
parameters may have converged before considering the highest possible order. After training 
the model, we can split these compressed parameters into the original ones as needed to 
make predictions for test cases. We show in detail how to compress parameters for logistic 
sequence prediction and logistic classification models. Experiments on both simulated and 
real data demonstrate that a huge number of parameters can indeed be reduced by our 
compression method. 

63 
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3.1 Introduction 

In many regression and classification problems, the response variable y depends on high- 
order interactions of "features" (also called "covariates" , "inputs" , "predictor variables" , or 
"explanatory variables"). Some complex human diseases are found to be related to high- 
order interactions of susceptibility genes and environmental exposures (Ritchie et. al. 2001). 
The prediction of the next character in English text is improved by using a large number 
of preceding characters (Bell, Cleary and Witten 1990). Many biological sequences have 
long-memory properties. 

When the features are discrete, we can employ high-order interactions in regression and 
classification models by introducing, as additional predictor variables, the indicators for each 
possible interaction pattern, equal to 1 if the pattern occurs for a subject and otherwise. 
In this chapter we will use "features" for the original discrete measurements and "predictor 
variables" for these derived variables, to distinguish them. The number of such predictor 
variables increases exponentially with the order of interactions. The total number of order-fc 
interaction patterns with k binary (0/1) features is 2 k , accordingly we will have 2 k predictor 
variables. A model with interactions of even a moderate order is prohibitive in real applica- 
tions, primarily for computational reasons. People are often forced to use a model with very 
small order, say only 1 or 2, which, however, may omit useful high-order predictor variables. 

Besides the computational considerations, regression and classification with a great many 
predictor variables may "overfit" the data. Unless the number of training cases is much larger 
than the number of predictor variables the model may fit the noise instead of the signal in 
the data, with the result that predictions for new test cases are poor. This problem can 
be solved by using Bayesian modeling with appropriate prior distributions. In a Bayesian 
model, we use a probability distribution over parameters to express our prior belief about 
which configurations of parameters may be appropriate. One such prior belief is that a 
parsimonious model can approximate the reality well. In particular, we may believe that most 
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high-order interactions are largely irrelevant to predicting the response. We express such a 
prior by assigning each regression coefficient a distribution with mode 0, such as a Gaussian 
or Cauchy distribution centered at 0. Due to its heavy tail, a Cauchy distribution may be 
more appropriate than a Gaussian distribution to express the prior belief that almost all 
coefficients of high order interactions are close to 0, with a very small number of exceptions. 
Additionally, the priors we use for the widths of Gaussian or Cauchy distributions for higher 
order interaction should favor small values. The resulting joint prior for all coefficients favors 
a model with most coefficients close to 0, that is, a model emphasizing low order interactions. 
By incorporating such prior information into our inference, we will not overfit the data with 
an unnecessarily complex model. 

However, the computational difficulty with a huge number of parameters is even more 
pronounced for a Bayesian approach than other approaches, if we have to use Markov chain 
Monte Carlo methods to sample from the posterior distribution, which is computationally 
burdensome even for a moderate number of parameters. With more parameters, a Markov 
chain sampler will take longer for each iteration and require more memory, and may need 
more iterations to converge or get trapped more easily in local modes. Applying Markov chain 
Monte Carlo methods to regression and classification with high-order interactions therefore 
seems infeasible. 

In this chapter, we show how these problems can be solved by effectively reducing the 
number of parameters in a Bayesian model with high-order interactions, using the fact that 
in a model that uses all interaction patterns, from a low order to a high order, many predictor 
variables have the same values for all training cases. For example, if an interaction pattern 
occurs in only one training case, all the interaction patterns of higher order contained in it 
will also occur in only that case and have the same values for all training cases — 1 for that 
training case and for all others. Consequently, only the sum of the coefficients associated 
with these predictor variables matters in the likelihood function. We can therefore use a 
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single "compressed" parameter to represent the sum of the regression coefficients for a group 
of predictor variables that have the same values in training cases. For models with very high 
order of interactions, the number of such compressed parameters will be much smaller than 
the number of original parameters. If the priors for the original parameters are symmetric 
stable distributions, such as Gaussian or Cauchy, we can easily find the prior distributions 
of these compressed parameters, as they are also symmetric stable distributions of the same 
type. In training the model with Markov chain Monte Carlo methods we need to deal only 
with these compressed parameters. After training the model, the compressed parameters 
can be split into the original ones as needed to make predictions for test cases. Using our 
method for compressing parameters, one can handle Bayesian regression and classification 
problems with very high order of interactions in a reasonable amount of time. 

This chapter will be organized as follows. We first describe Bayesian logistic sequence pre- 
diction models and Bayesian logistic classification models to which our compression method 
can be applied. Then, in Section 3.3 we describe in general terms the method of compressing 
parameters, and how to split them to make predictions for test cases. We then apply the 
method to logistic sequence models in Section 3.4, and to logistic classification models in 
Section 3.5. There, we will describe the specific schemes for compressing parameters for 
these models, and use simulated data and real data to demonstrate our method. We draw 
conclusions and discuss future work in Section 3.6. 

3.2 Two Models with High-order Interactions 

3.2.1 Bayesian Logistic Sequence Prediction Models 

We often need to predict the next state of a sequence given its preceding states, for example 
in speech recognition (Jelinek 1998), in text compression (Bell, Cleary, and Witten 1990), 
and in many others. We write a sequence of length O + 1 as x±, . . . ,xq,xq+i, where x t 
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Figure 3.1: A picture of the coefficients, /3, for all patterns in binary sequences of length 
= 3. /3[A 1 A 2 A- i ] is associated with the pattern written as [A1A2-A3], with A t — meaning 
that x t is allowed to be either 1 or 2, in other words, x t is ignored in defining this pattern. 
For example, /3[ oo] is the intercept term. These coefficients are used in defining the linear 
function I ((xi, x 2 , £3), (3) in the logistic model (3.1). For each combination of (xi,x 2 ,x 3 ) on 
the left column, / ((xi, x 2 , £3), /3) is equal to the sum of /3's along the path linked by lines, 
from P[ XlX2Xa ] to /3[ooo]- 



takes values from 1 to K t , for t = 1, . . . , O, and xo+i takes values from 1 to K. We call 
x\,...,xo = xi-o the historic sequence. For subject i we write its historic sequence and 
response as x^) Q and Iq+i- We are interested in modelling the conditional distribution 
P(x 0+ i I x 1:0 ). 

An interaction pattern V is written as [^y^ . . . Aq], where A t can be from to K t , with 
At = meaning that Xt can be any value from 1 to K t . For example, [0 . . . 01] denotes the 
pattern that fixes xo — 1 and allows x±, . . . ,xo-i to be any values in their ranges. When 
all nonzero elements of V are equal to the corresponding elements of a historic sequence, 
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X\ : o, we say that pattern V occurs in Xi-o, or pattern V is expressed by X\ : o, denoted by 
X\.o G V. We will use the indicator I{x\,o G V) as a predictor variable, whose coefficient 
is denoted by /3p. For example, /3[o...o] is the intercept term. A logistic model assigns each 

(k) 

possible value of the response a linear function of the predictor variables. We use j3 v to 
denote the coefficient associated with pattern V and used in the linear function for xo+i = k. 

For modeling sequences, we consider only the patterns where all zeros (if any) are at the 
start. Let us denote all such patterns by «S. We write all coefficients for xo+i = k, i.e., 
(3^ | V G <sj, collectively as /3^ k \ Figure (3.1) displays for binary sequence of length 
= 3, for some k, placed in a tree-shape. 

Conditional on f3^\ . . . , (3^ and Xi : q, the distribution of xq+i is defined as 



exp(/ (aWg)) 
Ef =1 ex P (/ (x 1:0 , 



P(x 0+ i = k | x l:Q , /3 (1) , • • • , ffi*>) = "^'ZT" 'Ls (3- 1 



where 



o 



l(x,. ,^) = E^ )j (^^)=^o] + EAoV,o] (3-2) 
veS 



t=i 



In Figure 3.1, we display the linear functions for each possible combination of (x\,X2, £3) 
on the left column, by linking together all /3's in the summation (3.2) with lines, from (3[ XlX2X3 ] 

to (3[ooo] ■ 

The prior for each (3^ is a Gaussian or Cauchy distribution centered at 0, whose width 
depends on the order, o(V), of V, which is the number of nonzero elements of V. There are 
0+1 such width parameters, denoted by (To, ... , oo- The cr D 's are treated as hyperparameters, 
assigned Inverse Gamma prior distributions with some shape and rate parameters, leaving 
their values to be determined by the data. In summary, the hierarchy of the priors is: 
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cr D ~ Inverse-Gamma (a Q , (a Q + 1) w ), for o — 0, . . . , O 

(k) (3 - 3) 
$P | cr o(P ) ~ Cauchy(0,cr o(P) ) or N(0,a^ v) ), for P G 5 

where Inverse-Gamma (a, A) denotes an Inverse Gamma distribution with density function 
x -ot-i ya exp(— A/x)/r(o;). We express a and A in (3.3) so that the mode of the prior is w Q . 



3.2.2 Remarks on the Sequence Prediction Models 

The Inverse Gamma distributions have heavy upward tails when a is small, and particularly 
when a < 1, they have infinite means. An Inverse Gamma distribution with a Q < 1 and 
small w Q , favors small values around w , but still allows o~ to be exceptionally large, as 
needed by the data. Similarly, the Cauchy distributions have heavy two-sided tails. The 
absolute value of a Cauchy random variable has infinite mean. When a Cauchy distribution 
with center and a small width is used as the prior for a group of parameters, such as all 
/5's of the interaction patterns with the same order in (3.3), a few parameters may be much 
larger in absolute value than others in this group. As the priors for the coefficients of high- 
order interaction patterns, the Cauchy distributions can therefore express more accurately 
than the Gaussian distributions the prior belief that most high-order interaction patterns 
are useless in predicting the response, but a small number may be important. 

It seems redundant to use a (3 {k) for each k = 1, . . . , K in (3.1) since only the differences 
between (3^ matter in (3.1). A non-Bayesian model could fix one of them, say /3 , all equal 
to 0, so as to make the parameters identifiable. However, when K ^ 2, forcing f3^ = in 
a Bayesian model will result in a prior that is not symmetric for all k, which we may not 
be able to justify. When K = 2, we do require that (3^ are all equal to 0, as there is no 
asymmetry problem. 

Inclusion of (3-p other than the highest order is also a redundancy, which facilitates the 
expression of appropriate prior beliefs. The prior distributions of linear functions of similar 
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historic sequences x\-o are positively correlated since they share some common /3's, for exam- 
ple, in the model displayed by Figure 3.1, I ((1, 1, l),/3) and I ((2, 1, share /3[on] , P[ooi] 
and /3[ooo]- Consequently, the predictive distributions of xo are similar given similar X\-.o- By 
incorporating such a prior belief into our inference, we borrow "statistical strength" for those 
historic sequences with few replications in the training cases from other similar sequences 
with more replications, avoiding making an unreasonably extreme conclusion due to a small 
number of replications. 

3.2.3 Bayesian Logistic Classification Models 

In this section we define the general Bayesian classification models with high-order interac- 
tions. We write the feature vector of dimension p as (xx, . . . , x p ), or collectively a? 1:p , where 
Xt takes values from 1 to K t , for t — 1, . . . ,p. In this thesis, we consider only classification 
problems in which the response y is discrete, assumed to take values from 1 to K. But 
our compression method could be applied to regression problems without any difficulty The 
features and response for subject i are written as x^ and y( % \ We are interested in modeling 
the conditional distribution P(y \ Xi :p ). 

A pattern V is written as . . . A p ], where A t can be from to K t , with At — 

meaning that x t can be any value from 1 to K t . For example, [0 ... 01] denotes the pattern 
that fixes x p = 1 and allows x\, . . . , x p -\ to be any values in their ranges. When all nonzero 
elements of V are equal to the corresponding elements of a feature vector, Xi- p , we say that 
pattern V occurs in Xi, p , or pattern V is expressed by X\ xp , denoted by X\ :p G V . The 
number of nonzero elements of a pattern V is called the order of V, denoted by o(V). All 
the patterns of order o are denoted by T* , and all the patterns from order to order O are 
denoted by V 0:O = Uo=o V ° ■ AU the patterns that are of order o and expressed by a feature 
vector x 1:p , are denoted by T >0 Xvp — { [A\, ■ ■ ■ , A p ] \ A t = or x t and Y^t=i H-^t 7^ 0) = o}. 
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Figure 3.2: A picture displaying all the interaction patterns, from order to order 3, of 3 
features x±,X2,X3, where x t is either 1 or 2. A pattern, written as [A1A2-A3], is shown by 3 
numbers linked by lines, where A t can be an integer from to 2, with A t = meaning x t 
could be either 1 or 2. The order of Ai,A 2 ,A 3 on the above graph can be changed to any 
permutation of Ai, A 2 , A 3 . 
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There are totally y P J patterns in P£ lV For example, "P(i,2,i) = {[0,2,1], [1,0,1], [1,2,0]}. 

Figure 3.2 displays P 013 for 3 binary (1 / 2) features xi,X2,x^. The patterns expressed by a 
feature vector (xi,X2,x^) can be found from such a graph, by searching from the root along 
the lines pointing to A t = and A t = x t . 

We will use the indicator I{x\- P G V) as a predictor variable, with coefficient denoted by 
f3-p. For example, (3\q...q] is the intercept term. A logistic model assigns each possible value of 
the response a linear function of the predictor variables. We use to denote the coefficient 
associated with pattern V and used in the linear function for y = k. All the coefficients for 
y = k are written as f3 {k) = {$ ] \ V G V°'°}. 

A Bayesian logistic classification model using all interaction patterns from order to 
order O is defined as follows: 



P(v = k \ . . .,/**>) = ^ggC3gq» (3.4) 



where 



o 



l(x^,0^)= ® ) I{*i*ZP)=Y, E (3-5) 



The priors for are given in the same way as in (3.3). 

The remarks regarding the Bayesian sequence prediction models in Section 3.2.2 still 
apply to the above classification models. Compared with the classification models, the 
Bayesian sequence prediction models are more restrictive models, using only the interaction 
patterns with all zeros at the start as predictor variables. 
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3.3 Our Method for Compressing Parameters 

In this section we describe in general terms our method for compressing parameters in 
Bayesian models, and how the original parameters can later be retrieved as needed for use 
in making predictions for test cases. 

3.3.1 Compressing Parameters 

In the above high-order models, the regression parameters of the likelihood function can 
be divided into a number of groups such that the likelihood function depends only on the 
sums over these groups, as shown by equation (3.8) below. We first use the Bayesian logistic 
sequence models to illustrate this fact. The likelihood function of (3^ k \ for k — 1, . . . , K, is 
the product of probabilities in (3.1) applied to the training for i = 1, . . . , N 

(collectively denoted by T>). It can be written as follows: 



(When K = 2, flS ' is fixed at 0, and therefore not included in the above likelihood function. 
But for simplicity, we do not write another expression for K = 2.) 

As can be seen in (3.2), the function I (x\ : o, (3) is the sum of the /3's associated with the 
interaction patterns expressed by X\-o- If a group of interaction patterns are expressed by 
the same training cases, the associated /3's will appear simultaneously in the same factors 
of (3.6). The likelihood function (3.6) therefore depends only on the sum of these /3's, rather 
than the individual ones. Suppose the number of such groups is G. The parameters in group 
g are rewritten as @ g i, . . . , fl gt n g , and the sum of them is denoted by s g : 



N 



exp(l(aff ,/3<'o+i>)) 
£f =1 exp(Z(a^ ,/3W)) 



lA/3« ...,0W|2>)=I] 



(3.6) 



n, 



u 




for g = 1, ... , 



G 



(3.7) 



fc=i 
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The likelihood function can then be rewritten as: 

£ (#11, • • • ,/?l,ni, ••• , Pgi, ■ ■ ■ , Pcnc) 

■■■> X>g* =£(*i, ••• ,*?) (3.8) 

fc=l k=l J 

(The above /3's are only the regression coefficients for the interaction patterns occurring in 
training cases. The predictive distribution for a test case may use extra regression coefficients, 
whose distributions depend only on the priors given relevant hyperparameters.) 

We need to define priors for the (5 gk in a way that lets us easily find the priors of the 
s g . For this purpose, we could assign each (3 gk a symmetric stable distribution centered at 
with width parameter a gk . Symmetric stable distributions (Feller 1966) have the following 
additive property: If random variables Xi,...,X n are independent and have symmetric 
stable distributions of index a, with location parameters and width parameters o"i, . . . , a n , 
then the sum of these random variables, 5^" =1 ^Q, also has a symmetric stable distribution 
of index a, with location parameter and width parameter (5^=i cr") 1 ■ Symmetric stable 
distributions exist and are unique for a G (0,2]. The symmetric stable distributions with 
a = 1 are Cauchy distributions. The density function of a Cauchy distribution with location 
parameter and width parameter a is [ira(l + x 2 /a 2 )}~ 1 . The symmetric stable distributions 
with a = 2 are Gaussian distributions, for which the width parameter is the standard 
deviation. Since the symmetric stable distributions with a other than 1 or 2 do not have 
closed form density functions, we will use only Gaussian or Cauchy priors. That is, each 
parameter f3 g k has a Gaussian or Cauchy distribution with location parameter and width 
parameter a gk : 

(3 gk ~N(0,o- 2 gk ) or p gk ~ Cauchy(0, a gk ) (3.9) 
As can been seen from the definitions of the priors (equation (3.3)), some a gk may be common 
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Figure 3.3: A picture depicting the sampling procedure after compressing parameters. 

for different (3 g k, but for simplicity we denote them individually. We might also treat the 
CgkS as unknown hyperparameters, but again we assume them fixed for the moment. 

If the prior distributions for the /3 g kS are as in (3.9), the prior distribution of s g can be 
found using the property of symmetric stable distributions: 

( n 9 \ / n 3 \ 

0, XX^J or * fl ~ Cauchy ( 0, X> fffc J (3-10) 

Let us denote the density of s g in (3.10) by Pg (either a Gaussian or Cauchy), and denote 
s\,...,sg collectively by s. The posterior distribution can be written as follows: 

P(s \ V) = L(s u ... ,s G )P{( Sl ) •■■ P°(s G ) (3.11) 

where T> is the training data, and c(T>) is the marginal probability or density function of T>. 

Since the likelihood function L(s\, . . . , sq) typically depends on s±, . . . , sq in a compli- 
cated way, we may have to use some Markov chain sampling method to sample for s from 
distribution (3.11). 

3.3.2 Splitting Compressed Parameters 

After we have obtained samples of s g , we may need to split them into their original compo- 
nents P g i, . . . , P g ,n g to make predictions for test cases. This "splitting" distribution depends 
only on the prior distributions, and is independent of the training data T>. In other words, the 
splitting distribution is just the conditional distribution of /? s i, . . . , (3 gng given Y^jkLi figk — s g> 
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whose density function is: 

P(Pgl, ■ ■ .,Pg,n g -l I Sg) = 

where P g k is the density function of the prior for (3 gk . Note that f3 g ^ ng is omitted since it is 
equal to s g - Yl=i Pgk- 

As will be discussed in the Section 3.3.4, sampling from (3.12) can be done efficiently by a 
direct sampling method, which does not involve costly evaluations of the likelihood function. 
We need to use Markov chain sampling methods and evaluate the likelihood function only 
when sampling for s. Figure 3.3 shows the sampling procedure after compressing parameters, 
where f3 is a collective representation of f3 g k, for g = 1, . . . , G, k = 1, . . . , n g — 1. When we 
consider high-order interactions, the number of groups, G, will be much smaller than the 
number of /3 g kS. This procedure is therefore much more efficient than applying Markov chain 
sampling methods to all the original f3 g h parameters. 

Furthermore, when making predictions for a particular test case, we actually do not 
need to sample from the distribution (3.12), of dimension n g — 1, but only from a derived 
1-dimensional distribution, which saves a huge amount of space. 

Before discussing how to sample from (3.12), we first phrase this compressing-splitting 
procedure more formally in the next section to show its correctness. 

3.3.3 Correctness of the Compressing-Splitting Procedure 

The above procedure of compressing and splitting parameters can be seen in terms of a 
transformation of the original parameters (3 g k to a new set of parameters containing s 5 's, as 
defined in (3.7), in light of the training data. The posterior distribution (3.11) of s and the 
splitting distribution (3.12) can be derived from the joint posterior distribution of the new 
parameters. 



n a -L 



Pgk(Pgk) 



fc=l 



P, 



(3.12) 



fc=i 
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The invertible mappings from the original parameters to the new parameters are shown 
as follows, for g — 1, . . . , G, 

n g 

(A?l, • • • , Pg,n B -l, Pg,n g ) =>" (Pgl, Pg,n g -1, ^ Pgk) = (Pgl, • • • , Pg,n g -1, S g ) (3.13) 

k=l 

In words, the first n g — 1 original parameters /3 g kS are mapped to themselves (we might use 
another set of symbols, for example b g k, to denote the new parameters, but here we still 
use the old ones for simplicity of presentation while making no confusion), and the sum of 
all /? 9i fc's, is mapped to s g . Let us denote the new parameters j3 g k, for g — 1, . . . , G, k — 
1, . . . , n g — 1, collectively by (3, and denote si, . . . , s g by s. (Note that (3 does not include 
/3 g ,n g , for 9 — 1, ■ ■ ■ ,G. Once we have obtained the samples of s and f3 we can use /3 9i „ 9 = 
s g ~ Y^kLi Pgk to obtain the samples of (3 g ,n g -) 

The posterior distribution of the original parameters, j3 g k, is: 



ii 



> • • • 5 rCr,n G 



PG,r. 



1 / m »G \ G n g 

p ) = e fa> e ^ n n p #w*) 

^ ' \k=\ k=l / g=l k=l 



(3.14) 



By applying the standard formula for the density function of transformed random variables, 
we can obtain from (3.14) the posterior distribution of the s and (3: 



P{s,(3 | V) 



1 



c(V) 



L(s 1 , 



sa)f[ 

9=1 



n -l 



n p M 



gk) 



k=l 



Pg,n g \ s g 



n g -l 

E 

k=l 



det(J)|(3.15) 



where the |det(J)| is absolute value of the determinant of the Jacobian matrix, J, of the 
mapping (3.13), which can be shown to be 1. 

Using the additive property of symmetric stable distributions, which is stated in sec- 
tion 3.3.1, we can analytically integrate out f3 in P(s,/3 \ T>), resulting in the marginal 
distribution P(s \ V): 
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P(s | V) 



P(s,f3 | V)d(3 



(3.16) 



1 



c(V) 

G 

n 

9=1 
1 

W) 



L (>!, . . . , S G ) ■ 
n a -l 



k=l 



Ug-l 

£ 

k=l 



d/3 g i 



L( Sl ,...,s G )P°( Sl ) ••• P G (s G ) 



df3 g>ng 43.17) 
(3.18) 



The conditional distribution of (3 given T> and s can then be obtained as follows: 



P((3 | s,V) = P(s,f3 | P)/P(a | V) 



G 

n 

9=1 



n„— 1 



[ ^(A 



fc=l 



n — 1 



(3.19) 
(3.20) 



fe=i 



From the above expression, it is clear that P{(3 \ s,V) is unrelated to V, i.e., P(f3 \ s,T>) = 
P(/3 | s), and is independent for different groups. Equation (3.12) gives this distribution 
only for one group g. 



3.3.4 Sampling from the Splitting Distribution 

In this section, we discuss how to sample from the splitting distribution (3.12) to make 
predictions for test cases after we have obtained samples of si, . . . , s G . 

If we sampled for all the f3 g kS, storing them would require a huge amount of space when 
the number of parameters in each group is huge. We therefore sample for (3 conditional 
on s\, . . . , s G only temporarily, for a particular test case. As can be seen in Section 3.2.1 
and 3.2.3, the predictive function needed to make prediction for a particular test case, for 
example the probability that a test case is in a certain class, depends only on the sums of 
subsets of j3 g kS in groups. After re-indexing the /3 g kS in each group such that the . . . , /3 9ttg 
are those needed by the test case, the variables needed for making a prediction for the test 
case are: 
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4 = for<7 = l,...,G, (3.21) 



k=l 

When t g = 0, s* = 0, and when t g = n g , s g = s g . The predictive function may also use some 
sums of extra regression coefficients associated with the interaction patterns that occur in 
this test case but not in training cases. Suppose the extra regression coefficients need to be 
divided into Z groups, as required by the form of the predictive function, which we denote 
by . . . , Pi jU *, • • • , Pz,iy • • • 5 Pz,n* ■ The variables needed for making prediction for the test 
cases are: 



fc=i 



P* zk , for*=l,...,Z (3.22) 



In terms of the above variables, the function needed to make a prediction for a test case 
can be written as 



tl t G 



a ••• ,^2Pak, ■■■ ,52Pzk\ = a ( s v ■■■ ,**o,*;, ••• ,»%) (3-23) 



fc=i fc=i fc=i fe=i 



Let us write s\, . . . , Sq collectively as s*, and write s*, . . . , s* z as s*. The integral required 
to make a prediction for this test case is 

G 

/ a(s 4 ,s*) P(s*) P(s | X>) [ P(sJ I s 9 ) ds'rfs*. (3.24) 

^ 9=1 

The integral over s* is done by MCMC. We also need to sample for s* from P(s*), 
which is the prior distribution of s* given some hyperparameters (from the current MCMC 
iteration) and can therefore be sampled easily. Finally, we need to sample from P(s* | s g ), 
which can be derived from (3.12), shown as follows: 

P(4 I s g ) = P«(4) Pf( Ss - st) /P s g {s g ) (3.25) 
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where Pg 1 ^ and Pg Z> are the priors (either Gaussian or Cauchy) of ^2\ 9 (3 g k and Yllg+i Pgki 
respectively. We can obtain (3.25) from (3.12) analogously as we obtained the density of 
s g , that is, by first mapping (3 and s to a set of new parameters containing s and s*, 
then integrating away other parameters, using the additive property of symmetric stable 
distributions. The distribution (3.25) splits s g into two components. 

When the priors for the f3 g kS are Gaussian distributions, the distribution (3.25) is also a 
Gaussian distribution, given as follows: 



>(2) 



si I s n ~ N I s £ 



S? + S^ 



2 ' ^1 



1 - 



(3.26) 



where S? = Eti ^ and Y? 2 = £* 

, 9 +i°"|fc- Since (3-26) is a Gaussian distribution, we can 
sample from it by standard methods. 

When we use Cauchy distributions as the priors for the Pg^s, the density function of (3.25) 



is: 



c sf + (4)2 s| + (4 - ^ 



(3.27) 



where Si = Ylk=i a gk> ^2 = Yltg+i °~9 fc ' an< ^ ^ * s ^ ne normalizing constant given below 
by (3.29). 

When s g = and Si = S 2 , the distribution (3.27) is a t-distribution with 3 degrees 
of freedom, mean and width Si/v3, from which we can sample by standard methods. 
Otherwise, the cumulative distribution function (CDF) of (3.27) can be shown to be: 



F(4; Sg ,Si,S 5 



1 

c 



,i ^2 



r log 



p (arctan ( 



y2 



S 

Si 



Ps (arctan ^ 



(3.2* 
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where 



7r ( E i + S 2) , 
C ~ ^ ^ T3. I 7^ , ^ \2\ ' l d ^ y J 



E 1 E 2 (s 2 + (E 1 + E 2 ) 

S 9 



s g + 2 (E 2 + Si) S 2 + (s 2 _ S 2 )2 > 

.2 fv2 



(3.30) 



j - (g ~ ^2 

Si S 4 + 2 (E? + E|) ^ + (£?-£ 2 



Po = t= — o, (3.31) 



1 s 2 + (£ 2 - Si) 
„ = 3 v 1 27 (3 32) 

S 2 S 4 + 2 (S 2 + S 2) ,2 + (S 2 _ S 2 )2 

When s g ^ 0, the derivation of (3.28) uses the equations below from (3.33) to (3.35) as 
follows, where p = (a 2 — c)/b,q = b + q,r = pc — a 2 q, and we assume 4c — b 2 > 0, 



1 1 1 f x + p x + q 



x 2 + a 2 x 2 + bx + c r V x 2 + a 2 x 2 + bx + c 



(3.33) 



W +p 1 . 2 2 . p fX\H 

-du = - log(x + a ) H — arctan I — 1 + — (3.34) 



u 2 + a 2 2 a \aJ 2 

u + q j 1 2 L Iq — b ( 2x + b \ 7r 



-(iit = - log(x + 6x + c) H — ^=^= arctan I : + — (3.35) 



DC 



n 2 + 6n + c 2 &v y V4c - 6 2 \yjAc-b 2 

When s 9 = 0, the derivation of (3.28) uses the following equations: 



11 1/1 1 



x 2 + a 2 x 2 + b 2 b 2 — a 2 V a; 2 + a 2 x 2 + 6 2 



(3.36) 



m 2 + a 2 aV la/ 2 



du = — ('arctan f— H — (3.37) 



Since we can compute the CDF of (3.27) with (3.28) explicitly, we can use the inversion 
method to sample from (3.27), with the inverse CDF computed by some numerical method. 
We chose the Illinois method (Thisted 1988, Page 171), which is robust and fairly fast. 

When sampling for s\, . . . , Sq temporarily for each test case is not desired, for example, 
when we need to make predictions for a huge number of test cases at a time, we can still 
apply the above method that splits a Gaussian or Cauchy random variable into two parts 
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n g — 1 times to split s g into n g parts. Our method for compressing parameters is still useful 
because sampling from the splitting distributions uses direct sampling methods, which are 
much more efficient than applying Markov chain sampling method to the original parameters. 
However, we will not save space if we take this approach of sampling for all /3's. 

3.4 Application to Sequence Prediction Models 

In this section, we show how to compress parameters of logistic sequence prediction models 
in which states of a sequence are discrete, as defined in Section 3.2.1. To demonstrate our 
method, we use a binary data set generated using a hidden Markov model, and a data set 
created from English text, in which each state has 3 possibilities (consonant, vowel, and 
others). These experiments show that our compression method produces a large reduction 
in the number of parameters needed for training the model, when the prediction for the next 
state of a sequence is based on a long preceding sequence, i.e., a high-order model. We also 
show that good predictions on test cases result from being able to use a high-order model. 

3.4.1 Grouping Parameters of Sequence Prediction Models 

In this section, we describe a scheme for dividing the /3's into a number of groups, based on 
the training data, such that the likelihood function depends only on the sums in groups, as 
shown by (3.8). 

Since the linear functions for different values of response have the same form except the 
superscript, the way we divide (3^ into groups is the same for all k. Our task is to find the 
groups of interaction patterns expressed by the same training cases. 

Let us use E-p to denote the "expression" of the pattern V — the indices of training 
cases in which V is expressed, a subset of 1, . . . , N. For example, Ejo-o] = {1, • • • , N}. In 
other words, the indicator for pattern V has value 1 for the training cases in E-p, and for 
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others. We can display E-p in a tree-shape, as we displayed (3-p. The upper part of Figure 3.4 
shows such expressions for each pattern of binary sequence of length = 3, based on 3 
training cases: x^l = (1, 2, l),x^l = (2,1,2) and x^l = (1,1,2). From Figure 3.4, we can 
see that the expression of a "stem" pattern is equal to the union of the expressions of its 
"leaf" patterns, for example, -Epoo] — -^[ooi] U -^[002] • 

When a stem pattern has only one leaf pattern with non-empty expression, the stem and 
leaf patterns have the same expression, and can therefore be grouped together. This grouping 
procedure will continue by taking the leaf pattern as the new stem pattern, until encountering 
a stem pattern that "splits" , i.e. has more than one leaf pattern with non-empty expression. 
For example, -Ejooi], -^[021] and £'[121] in Figure 3.4 can be grouped together. All such patterns 
must be linked by lines, and can be represented collectively with a "superpattern" SP, 
written as [0 • • • 0A b ■ ■ ■ A ]f = {J b t=f [0---0Af-A o }, where 1 < b < f < O + 1, and 
in particular when t = + 1, [0 • • • 0A t ■ ■ ■ Ao] = [0 ■ ■ • 0]. One can easily translate the 
above discussion into a computer algorithm. Figure 3.5 describes the algorithm for grouping 
parameters of Bayesian logistic sequence prediction models, in a C-like language, using a 
recursive function. 

An important property of our method for compressing parameters of sequence prediction 
models is that given N sequences as training data, conceivably of infinite length, denoted by 
. . . , for i — 1, . . . , N, the number of superpatterns with unique expressions, and 
accordingly the number of compressed parameters, will converge to a finite number as O 
increases. The justification of this claim is that if we keep splitting the expressions following 
the tree shown in Figure 3.4, at a certain time, say t, every expression will be an expression 
with only 1 element (suppose we in advance remove the sequences that are identical with 
another one). When considering further smaller t, no more new superpattern with different 
expressions will be introduced, and the number of superpatterns will not grow. The number 
of the compressed parameters, the regression coefficients for the superpatterns, will therefore 
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Figure 3.4: A picture showing that the interaction patterns in logistic sequence prediction 
models can be grouped, illustrated with binary sequences of length = 3, based on 3 training 
cases shown in the upper-right box. E-p is the expression of the pattern (or superpattern) V 
- the indices of the training cases in which the V is expressed, with meaning the empty 
set. We group the patterns with the same expression together, re-represented collectively 
by a "superpattern", written as [0 ■ ■ -0A b ■ ■ ■ A Q ]f, meaning |jf =6 [0 • • -0A t ■ ■ -Ao], where 
1 < b < f < O + 1, and in particular when t = O + 1, [0 • • • 0A t ■ • ■ A ] = [0 • • • 0]. We also 
remove the patterns not expressed by any of the training cases. Only 5 superpatterns with 
unique expressions are left in the lower picture. 
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INPUTS: 

N: number of training cases 

O: length of sequences 
X: training data, N X O matrix 
OUTPUTS: 

LSP: list of superpatterns 
LE: list of expressions 

INPUTS of 'DIVERGE' (shown on the right): 

E: expression, subset of 1, ... ,N 

SP: superpattern, structure with members 
P: pattern, vector of length O 

f : index of fixed state in P 

ALGORITHM: 
E= {1, ... ,N} 

SP.f = 0+ 1, SP.P=(0, ... ,0) 
LSP = NULL 
LE = NULL 
DIVERGE(E,SP) 
RETURN LE and LSP 



DIVERGE(E, SP) 
{ 

b = SP.f- 1 

while ( # of unique values in X[E][b] is 1 & b >0) 
{ SP.P [b] = the unique value in X[E] [b] 
b = b- 1 

1 

Add SP to LSP 
Add E to LE 
if ( b > ) 

{ for( x in unique values in X[E][b] ) 
{ SubE={iinE:X[i][b]=x} 
Set NSP = SP 
NSP.f = b , NSP.P[b] = x 
DIVERGE(SubE, NSP) 

} 

} 



Figure 3.5: The algorithm for grouping parameters of Bayesian logistic sequence prediction 
models. To group parameters, we call function "DIVERGE" with the initial values of ex- 
pression E = {1,...,N} and superpattern SP = [0...0]o+i, as shown in above picture, 
resulting in two lists of the same length, LE and LSP, respectively storing the expressions 
and the corresponding superpatterns. Note that the first index of an array is assumed to be 
1 , and that the X [E] [b] means a 1-dimension subarray of X in which the row indices are in 
E and the column index equals b. 
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not grow after the time t. 

In contrast, after the time t when each interaction pattern is expressed by only 1 training 
case, if the order is increased by 1, the number of interaction patterns is increased by the 
number of training cases. The regression coefficients associated with these original inter- 
action patterns, called the original parameters thereafter, will grow linearly with the order 
considered. Note that these original parameters do not include the regression coefficients for 
those interaction patterns not expressed by any training case. The total number of regression 
coefficients defined by the model grows exponentially with the order considered. 

3.4.2 Making Prediction for a Test Case 

Given . . . , (3^ K \ the predictive probability for the next state x* Q+1 of a test case for 
which we know the historic sequence x\. Q can be computed using equation (3.1), applied to 
x\. . A Monte Carlo estimate of P(x* 0+1 = k | x\. ,T>) can be obtained by averaging (3.1) 
over the Markov chain samples from the posterior distribution of j3^\ . . . , j3^ K \ 

Each of the + 1 patterns expressed by the test is either expressed by some 

training case (and therefore belongs to one of the superpatterns), or is a new pattern (not 
expressed by any training case). Suppose we have found 7 superpatterns. The + 1 /3's in 
the linear function l(x\. , (3^) can accordingly be divided into 7 + 1 groups (some groups 
may be empty). The function Z(x*. , /?( fc )) can be written as the sum of the sums of the 
j3 7 s over these 7 + 1 groups. Consequently, P(x* 0+1 = k | x\. Q ) can be written in the form 
of (3.23). As discussed in Section 3.3.4, we need to only split the sum of the /3's associated 
with a superpattern, i.e., a compressed parameter s g , into two parts, such that one of them 
is the sum of those (3 expressed by the test using the splitting distribution (3.25). 

It is easy to identify the patterns that are also expressed by x\. Q from a superpattern 
[0---A b ---A o ]f. If( x *fi ■ ■ ■ 1 x *o) 7^ (Af > • • • ' A-o)-, none of the patterns in [0 • • • • • • Aq]/ 

are expressed by x\. , otherwise, if . . . , x* Q ) = {Ay, . . . , Aq) for some b' {b < b' < f) , all 
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patterns in [0 • • ■ Ay ■ ■ ■ Ao]f are expressed by x\. . 

3.4.3 Experiments with a Hidden Markov Model 

In this section we apply Bayesian logistic sequence prediction modeling, with or without our 
compression method, to data sets generated using a Hidden Markov model, to demonstrate 
our method for compressing parameters. The experiments show that when the considered 
length of the sequence O is increased, the number of compressed parameters will converge 
to a fixed number, whereas the number of original parameters will increase linearly. Our 
compression method also improves the quality of Markov chain sampling in terms of auto- 
correlation. We therefore obtain good predictive performances in a small amount of time 
using long historic sequences. 

The Hidden Markov Model Used to Generate the Data 

Hidden markov models (HMM) are applied widely in many areas, for example, speech recog- 
nition (Baker 1975), image analysis (Romberg et.al. 2001), computational biology (Sun 
2006). In a simple hidden Markov model, the observable sequence {x t | t — 1, 2, . . .} is mod- 
eled as a noisy representation of a hidden sequence {h t | t — 1, 2, . . .} that has the Markov 
property (the distribution of ht given h t -i is independent with the previous states before 
/i 4 „i). Figure 3.6 displays the hidden Markov model used to generate our data sets, showing 
the transitions of three successive states. The hidden sequence h t is an Markov chain with 
state space {1, . . . , 8}, whose dominating transition probabilities are shown by the arrows in 
Figure 3.6, each of which is 0.95. However, the hidden Markov chain can move from any state 
to any other state as well, with some small probabilities. If h t is an even number, xt will be 
equal to 1 with probability 0.95 and 2 with probability 0.05, otherwise, x t will be equal to 2 
with probability 0.95 and 1 with probability 0.05. The sequence {xt | t = 1, 2, . . .} generated 
by this exhibits high-order dependency, though the hidden sequence is only a Markov chain. 
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We can see this by looking at the transitions of observable Xt in Figure 3.6. For example, if 
x\ — 1 (rectangle) and X2 = 2 (oval), it is most likely to be generated by hi = 2 and l%2 = 3, 
since this is the only strong connection from the rectangle to the oval, consequently, = 8 
is most likely to to be the next, and X3 is therefore most likely to be 1 (rectangle). 



hj h 2 h 3 




Figure 3.6: A picture showing a Hidden Markov Model, which is used to generate sequences 
to demonstrate Bayesian logistic sequence prediction models. Only the dominating transition 
probabilities of 0.95 are shown using arrows in the above graph, while from any state the 
hidden Markov chain can also move to any other state with a small probability. When ht is 
in a rectangle, Xt is equal to 1 with probability 0.95, and 2 with probability 0.05, otherwise, 
when h t is in an oval, x t is equal to 2 with probability 0.95, and 1 with probability 0.05. 

3.4.4 Specifications of the Priors and Computation Methods 
The Priors for the Hyperprameters 

We fix o"o at 5 for the Cauchy models and 10 for the Gaussian models. For o > 0, the prior 
for a is Inverse Gamma(a D , [a Q + l)w ), where a Q and w Q are: 
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a o = 0.25, w = 0.1/o, foro=l,...,0 (3.38) 
The quantiles of Inverse-Gamma(0.25, 1.25 x 0.1), the prior of a\, are shown as follows: 
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The quantiles of other a Q can be obtained by multiplying those of o\ by 1/ 0. 
The Markov Chain Sampling Method 

We use Gibbs sampling to sample for both the s 9 's (or the {3 g kS when not applying our 
compression method) and the hyperparameters, a . These 1-dimensional conditional distri- 
butions are sampled using the slice sampling method (Neal 2003), summarized as follows. 
In order to sample from a 1-dimensional distribution with density f(x), we can draw points 
(x, y) from the uniform distribution over the set {(x, y) | < y < f(x)}, i.e., the region of 
the 2-dimensional plane between the x-axis and the curve of f(x). One can show that the 
marginal distribution of x drawn this way is f{x). We can use Gibbs sampling scheme to 
sample from the uniform distribution over {(x, y) | < y < f{x)}. Given x, we can draw y 
from the uniform distribution over {y | < y < f(x)}. Given y, we need to draw x from the 
uniform distribution over the "slice", S = {x \ f(x) > y}. However, it is generally infeasible 
to draw a point directly from the uniform distribution over S. (Neal 2003) devises several 
Markov chain sampling schemes that leave this uniform distribution over 5* invariant. One 
can show that this updating of x along with the previous updating of y leaves f(x) invariant. 
Particularly we chose the "stepping out" plus "shrinkage" procedures. The "stepping out" 
scheme first steps out from the point in the previous iteration, say Xo, which is in S, by 
expanding an initial interval, /, of size w around xq on both sides with intervals of size w, 
until the ends of / are outside S, or the number of steps has reached a pre-specified number, 
m. To guarantee correctness, the initial interval I is positioned randomly around xq, and 



3 Compressing Parameters in Bayesian Models with High- order Interactions 



90 



m is randomly aportioned for the times of stepping right and stepping left. We then keep 
drawing a point uniformly from / until obtaining an x in S. To facilitate the process of 
obtaining an x in S, we shrink the interval / if we obtain an x not in S by cutting off the 
left part or right part of I depending on whether x < Xq or x > Xq. 

We set w = 20 when sampling for /3's if we use Cauchy priors, considering that there might 
be two modes in this case, and set w = 10 if we use Gaussian priors. We set w — 1 when 
sampling for a Q . The value of m is 50 for all cases. We trained the Bayesian logistic sequence 
model, with the compressed or the original parameters, by running the Markov chain 2000 
iterations, each updating the /3's 1 time, and updating the cr's 10 times, both using slice 
sampling. The first 750 iterations were discarded, and every 5th iteration afterward was 
used to predict for the test cases. The number of 750 was chosen empirically after looking 
at many trial runs of Markov chains for many different circumstances. 

The above specification of Markov chain sampling and the priors for the hyperparameters 
will be used for all experiments in this chapter, including the experiments on classification 
models discussed in Section 3.5. 

Experiment Results 

We used the HMM in Figure 3.6 to generate 5500 sequences with length 21. We used 5000 
sequences as test cases, and the remaining 500 as the training cases. We tested the prediction 
methods by predicting X21 based on varying numbers of preceding states, O, chosen from 
the set {1, 2, 3, 4, 5, 7, 12, 15, 17, 20}. 

Figure 3.7 compares the number of parameters and the times used to train the model, with 
and without our compression method. It is clear that our method for compressing parameters 
reduces greatly the number of parameters. The ratio of the number of compressed parameters 
to the number of the original ones decreases with the number of preceding states, O. For 
example, the ratio reaches 0.207 when O = 20. This ratio will reduce to when considering 
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even bigger O, since the number of original parameters will grow with O while the number 
of compressed parameters will converge to a finite number, as discussed in Section 3.4.1. 
There are similar reductions for the training times with our compression method. But the 
training time with compressed parameters will not converge to a finite amount, since the time 
used to update the hyperparameters (er 's) grows with order, O. Figure 3.7 also shows the 
prediction times for 5000 training cases. The small prediction times show that the methods 
for splitting Gaussian and Cauchy variables are very fast. The prediction times grow with 
O because the time used to identify the patterns in a superpattern expressed by a test case 
grows with O. The prediction times with the original parameters are not shown in Figure 3.7, 
since we do not claim that our compression method saves prediction time. (If we used the 
time-optimal programming method for each method, the prediction times with compressed 
parameters should be more than without compressing parameters since the method with 
compression should include times for identifying the patterns from the superpattern for test 
cases. With our software, however, prediction times with compression are less than without 
compression, which is not shown in Figure 3.7, because the method without compression 
needs to repeatedly read a huge number of the original parameters into memory from disk.) 

Compressing parameters also improves the quality of Markov chain sampling. Figure 3.8 
shows the autocorrelation plots of the hyperparameters a Q , for o = 10, 12, 15, 17, 20, when 
the length of the preceding sequence, O, is 20. It is clear that the autocorrelation decreases 
more rapidly with lag when we compress the parameters. This results from the compressed 
parameters capturing the important directions of the likelihood function (i.e. the directions 
where a small change can result in large a change of the likelihood). We did not take the time 
reduction from compressing parameters into consideration in this comparison. If we rescaled 
the lags in the autocorrelation plots according to the computation time, the reduction of 
autocorrelation of Markov chains with the compressed parameters would be much more 
pronounced. 
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Figure 3.7: Plots showing the reductions of the number of parameters and the training time 
with our compression method using the experiments on a data set generated by a HMM. The 
upper-left plot shows the number of the compressed and the original parameters based on 
500 training sequences for O — 1,2, 3, 4, 5, 7, 10, 12, 15, 17, 20, their ratios are shown in the 
upper-right plot. In the above plots, the lines with o are for the methods with parameters 
compressed, the lines with x are for the methods without parameters compressed, the dashed 
lines are for the methods with Gaussian priors, and the dotted lines are for the methods with 
Cauchy priors. The lower-left plot shows the training times for the methods with and without 
parameters compressed. The lower-right plot shows the prediction time only for the methods 
with parameters compressed. 
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Figure 3.8: The autocorrelation plots of cr 's for the experiments on a data set generated by 
a HMM, when the length of the preceding sequence O = 20. We show the autocorrelations 
of o~ , for o = 10, 12, 15, 17, 20. In the above plots, "Gaussian" in the titles indicates the 
methods with Gaussian priors, "Cauchy" indicates with Cauchy priors, "comp" indicates 
with parameters compressed, "no comp" indicates without parameters compressed. 
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Figure 3.9: Plots showing the predictive performance using the experiments on a data set 
generated by a HMM. The left plots show the error rates and the right plots show the 
average minus log probabilities of the true responses in the test cases. The upper plots 
show the results when using the Cauchy priors and the lower plots shows the results when 
using the Gaussian priors. In all plots, the lines with o are for the methods with parameters 
compressed, the lines with x are for the methods without parameters compressed. The 
numbers of the training and test cases are respectively 500 and 5000. The number of classes 
of the response is 2. 



3 Compressing Parameters in Bayesian Models with High- order Interactions 



95 



Finally, we evaluated the predictive performance in terms of error rate (the fraction of 
wrong predictions in test cases), and the average minus log probability (AMLP) of observing 
the true response in a test case based on the predictive probability for different classes. The 
performance of with and without compressing parameters are the same, as should be the 
case in theory, and will be in practice when the Markov chains for the two methods converge 
to the same modes. Performance of methods with Cauchy and Gaussian priors is also similar 
for this example. The predictive performance is improved when O goes from 1 to 5. When 
O > 5 the predictions are slightly worse than with O = 5 in terms of AMLP. The error 
rates for O > 5 are almost the same as for = 5. This shows that the Bayesian models 
can perform reasonably well even when we consider a very high order, as they avoid the 
overfitting problem in using complex models. We therefore do not need to restrict the order 
of the Bayesian sequence prediction models to a very small number, especially after applying 
our method for compressing parameters. 

3.4.5 Experiments with English Text 

We also tested our method using a data set created from an online article from the website 
of the Department of Statistics, University of Toronto. In creating the data set, we encoded 
each character as 1 for vowel letters (a,e,i,o,u), 2 for consonant letters, and 3 for all other 
characters, such as space, numbers, special symbols, and we then collapsed multiple occur- 
rences of "3" into only 1 occurrence. The length of the whole sequence is 3930. Using it we 
created a data set with 3910 overlaped sequences of length 21, and used the first 1000 as 
training data. 

The experiments were similar to those in Section 3.4.3, with the same priors and the 
same computational specifications for Markov chain sampling. Figures 3.10, 3.11, 3.12, and 
3.13 show the results. All the conclusions drawn from the experiments in Section 3.4.3 are 
confirmed in this example, with some differences in details. In summary, our compression 
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Figure 3.10: Plots showing the reductions of the number of parameters and the training and 
prediction time with our compression method using the experiments on English text. The 
upper-left plot shows the number of the compressed and the original parameters based on 
500 training sequences for O — 1,2, 3, 4, 5, 7, 10, 12, 15, 17, 20, their ratios are shown in the 
upper-right plot. In the above plots, the lines with o are for the methods with parameters 
compressed, the lines with x are for the methods without parameters compressed, the dashed 
lines are for the methods with Gaussian priors, and the dotted lines are for the methods with 
Cauchy priors. The lower-left plot shows the training times for the methods with and without 
parameters compressed. The lower-right plot shows the prediction time only for the methods 
with parameters compressed. 
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Figure 3.11: The autocorrelation plots of the cr 's for the experiments on English text data, 
when the length of the preceding sequence O = 20. We show the autocorrelation plot 
of <r , for o = 10, 12, 15, 17, 20. In the above plots, "Gaussian" in the titles indicates the 
methods with Gaussian priors, "Cauchy" indicates with Cauchy priors, "comp" indicates 
with parameters compressed, "no comp" indicates without parameters compressed. 
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Figure 3.12: Plots showing the predictive performance using the experiments on English 
text data. The left plots show the error rate and the right plots show the average minus log 
probability of the true response in a test case. The upper plots show the results when using 
the Cauchy priors and the lower plots shows the results when using the Gaussian priors. In 
all plots, the lines with o are for the methods with parameters compressed, the lines with x 
are for the methods without parameters compressed. The numbers of the training and test 
cases are respectively 1000 and 2910. The number of classes of the response is 3. 
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Figure 3.13: Scatterplots of medians of all compressed parameters, s, of Markov chain sam- 
ples in the last 1250 iterations, for the models with Cauchy and Gaussian priors, fitted with 
English text data, with the length of preceding sequence O = 10, and with the parameters 
compressed. The right plot shows in a larger scale the rectangle (—2, 2) x (—2, 2). 

method reduces greatly the number of parameters, and therefore shortens the training process 
greatly. The quality of Markov chain sampling is improved by compressing parameters. 
Prediction is very fast using our splitting methods. The predictions on the test cases are 
improved by considering higher order interactions. From Figure 3.12, at least some order 10 
interactions are useful in predicting the next character. 

In this example we also see that when Cauchy priors are used Markov chain sampling with 
the original parameters may have been trapped in a local mode, resulting in slightly worse 
predictions on test cases than with the compressed parameters, even though the models used 
are identical. 

We also see that the models with Cauchy priors result in better predictions than those 
with Gaussian priors for this data set, as seen from the plots of error rates and AMLPs. To 
investigate the difference of using Gaussian and Cauchy priors, we first plotted the medians 
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Figure 3.14: Plots of Markov chain traces of three compressed parameters (each contains only 
one (3) from experiments on English text with 10 preceding states, with Cauchy or Gaussian 
priors. In each plot, three different lines show three indepedent runs. The parameters are 
annotated by their original meanings in English sequence. For example, £ _CC:V stands 
for the parameter for predicting that the next character is a "vowel" given preceding three 
characters are "space, consonant, consonant". 
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of Markov chains samples (in the last 1250 iterations) of all compressed parameters, s, for 
the model with O = 10, shown in Figure 3.13, where the right plot shows in a larger scale 
the rectangle (—2,2) x (—2,2). This figure shows that a few (3 with large medians in the 
Cauchy model have very small corresponding medians in the Gaussian model. 

We also looked at the traces of some compressed parameters, as shown in Figure 3.14. 
The three compressed parameters shown all contain only a single (3. The plots on the top 
are for the (3 for "CC: V" , used for predicting whether the next character is a vowel given the 
preceding two characters are consonants; the plots in the middle are for "__CC:V", where 
denotes a space or special symbol; the plots on the bottom are for "CCVCVCC:V", 
which had the largest median among all compressed parameters in the Cauchy model, as 
shown by Figure 3.13. The regression coefficient (3 for "CC:V" should be close to by our 
common sense, since two consonants can be followed by any of three types of characters. 
We can very commonly see "CCV", such as "the", and "CC__", such as "with__", and not 
uncommonly see "CCC", such as "technique" , "world" , etc. The Markov chain trace of this j3 
with a Cauchy prior moves in a smaller region around than with a Gaussian prior. But if we 
look back one more character, things are different. The regression coefficient (3 for "__CC:V" 
is fairly large, which is not surprising. The two consonants in "__CC:V" stand for two letters 
in the beginning of a word. We rarely see a word starting with three consonants or a word 
consisting of only two consonants. The posterior distribution of this 9 for both Cauchy and 
Gaussian models favor positive values, but the Markov chain trace for the Cauchy model 
can move to much larger values than for the Gaussian model. As for the high-order pattern 
"CCVCVCC" , it matches words like "statistics" or "statistical" , which repeatedly appear in 
an article introducing a statistics department. Again, the Markov chain trace of this (3 for 
the Cauchy model can move to much larger values than for Gaussian model, but sometimes 
it is close to 0, indicating that there might be two modes for its posterior distribution. 

The above investigation reveals that a Cauchy model allows some useful (3 to be much 
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larger in absolute value than others while keeping the useless (3 in a smaller region around 
than a Gaussian model. In other words, Cauchy models are more powerful in finding the 
information from the many possible high-order interactions than Gaussian models, due to 
the heavy two-sided tails of Cauchy distributions. 



3.5 Application to Logistic Classification Models 

3.5.1 Grouping Parameters of Classification Models 

As we have seen in sequence prediction models, the regression coefficients for the patterns 
that are expressed by the same training cases can be compressed into a single parameter. 
We present an algorithm to find such groups of patterns. Our algorithm may not the only 
one possible and may not be the best. 

Our algorithm uses a "superpattern", SP, to represent a set of patterns with some 

common property, written as ( ^ ' ' ' ^ p , where A t is the pattern value for position t 

V h...I p h 

(an integer from to Kt), It is binary (0/1) with 1 indicating A t is fixed for this superpattern, 
/ is the number of fixed positions (ie / = Y^t=i ^t)> an d indicates the smallest order of all 
patterns in this superpattern, equal to the sum of nonzero values of those A t with I t = 1 
(i.e. o = J2t=i 1(1* = 7^ 0))- Such a superpattern represents the union of all patterns 
with order not greater than O, with values at the fixed positions (with I t = 1) being At, 
but the pattern values at unfixed positions (with I t = 0) being either or A t . For example, 

if O = 3, the superpattern ^ |qq^q ) ^ s com P ose d °^ ( g ) ' ( ^ ) ' an< ^ ( 2 ) P a ^ erns 
respectively of order 1, 2 and 3, listed as follows: 



order 1 : [10000] 

order 2 : [12000], [10300], [10004] 
order 3 : [10304], [12004], [12300] 



(3.39) 
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The algorithm is inspired by the display of all interaction patterns in Figure 3.2. It starts 



features (1, . . . ,p). After choosing a feature x t from the unconsidered features by the way 
described below, the current expression is split for each value of the x t , as done by the 
algorithm for sequence prediction models, additionally the whole expression is also passed 
down for the pattern with A t = 0. When we see that we can not split the expression by 
any of the unconsidered features, all the following patterns can be grouped together and 
represented with a superpattern. In Figure 3.15, we use training data with 3 binary features 
and only 3 cases to illustrate the algorithm. We give the algorithm in a C-like language in 
Figure 3.16, which uses a recursive function. 

In choosing a feature for splitting a expression, we look at the diversity of the values of 
the unconsidered features restricted on the current expression. By this way the expression 
is split into more expressions but each may be smaller. We therefore more rapidly reach 
the expression that can not be split further. The diversity of a feature x t restricted on the 
current expression is measured with the entropy of the relative frequency of the values of x t , 
i.e., — ^2iPi logtjOj), where p^s are the relative frequency of the possible values of the feature 
restricted on the expression. When two features have the same entropy value, we choose the 
one with smaller index t. Note that the entropy is always positive unless all the values of 
the feature x t restricted on the expression are the same. The resulting expressions found by 
this algorithm are not unique for each superpattern. 

In training the model, we need to compute the width of a parameter associated with a 
superpattern given the values of the hyperparameters cr 's. For Cauchy models, the width is 
equal to the sum of the hyperparameters of all patterns in the superpattern. For Gaussian 
models, the width is the square roots of the sum of the squares of the hyperparameters of all 
patterns in the superpattern. We therefore need only know the number of the patterns in 
the superpattern belonging to each order from to O. For a superpattern I 1 ' ' ' p ) , 



with the expression {1,...,N} for the superpattern 



( 



oo. ..o y 



and the unconsidered 
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Figure 3.15: This picture illustrates the algorithm for grouping the patterns of classification 
models using a training data of 3 binary features and 3 cases shown on the left-top corner. 
Starting from the expression for the intercept and all features in the unconsidered features 
set, we recursively split the current expression by the values of the feature with the biggest 
entropy (diversity) among the remaining unconsidered features. When the values of all 
remaining unconsidered features are the same for all training cases, for example, when the 
expression contains only one training case, all the following patterns can be grouped together. 
After grouping, all the expressions in dashed box are removed and grouped into their parent 
expressions, with the group of patterns being represented using a superpattern. 
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INPUTS: 

N: number of training cases 
p: number of features 
O: the order of the model (<=p) 
X: training data, N Xp matrix 

OUTPUTS: 

LSP: list of superpatterns 
LE: list of expressions 

INPUTS of 'DIVERGE' (shown on the right): 
E: expression, subset of 1, ... ,N 
SP: superpattern, structure with members 

P: pattern, array of 2 X p 

f : number of fixed positions 

o: smallest order 

FT: indice of unconsidered features 

ALGORITHM: 

E={1, ...,N} 

SP:o = 0,f=p,P=(0 0-0) 

FT = (l,2,...,p) 

LSP = NULL 
LE = NULL 
DIVERGE(E,SP,FT) 
RETURN LSP and LE 



DIVERGE(E, SP,FT) 

{ find the entropies of each column of X[E][FT] 
M = the biggest entropy 
mFT = index of the feature with entropy=M 

if(M == 0) 

{ Create new NSP, Set NSP= SP 

NSP.P[FT] = X[E[1]][FT], NSP.f -= # of FT 
Add NSP to LSP, Add E to LE 



else 



SubFT = FT with mFT removed 

if(# of SubFT == 0) 

{Add SP to LSP, Add E to LE} 

else DIVERGE(E,SP,SubFT) 

for(x in unique values of X[E][mFT]) 

{ Set NSP=SP 

NSP.P[mFT] = x, NSP.o ++ 
SubE = {i in E I X[i][mFT] =x} 

if(# of SubFT ==011 NSP.o == O) 

{Add SP to LSP, Add SubE to LE) 
else DIVERGE(subE,NSP,subFT) 



Figure 3.16: The algorithm for grouping the patterns of Bayesian logistic classification mod- 
els. To do grouping, we call the function "DIVERGE" with the initial values of expression 
E, superpattern SP and the unconsidered features FT shown as above, resulting in two lists 
of the same length, LE and LSP, respectively storing the expressions and the corresponding 
superpatterns. Note that the first index of an array are assumed to be 1, and thatXf-E] [FT] 
means sub matrix of X by restricting the rows in E and columns in FT. 3) The resulting 
expressions are not unique for each superpattern in LSP. We still need to merge those 
superpatterns with the same expression by directly comparing the expressions. 
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they are given as follows: 



# (patterns of order o + d) = < 




for d — 0, . . . , min(0 — o,p — f) 



(3.40) 



0. 



Otherwise 



V 



In predicting for a test case x*, we need to identify the patterns in a superpattern 



none of the patterns in the superpatterns are expressed by x*. Otherwise, if x* t ^ A t for any 
unfixed positions (with I t = 0) then the A t is set to 0. This results in a smaller superpattern, 
from which we can count the number of patterns belonging to each order using the formula 



For a fixed very large number of features, p, and a fixed number of cases, N, the number 
of the compressed parameters may have converged before considering the highest possible 
order, p. This can be verified by regarding the interaction patterns of a classification model as 
the union (non-exclusive) of the interaction patterns of p\ sequence models from permutating 
the indice of features. As shown for sequence models in Section 3.4.1, there is a certain order 
Oj, for j = 1, . . . ,p\, for each of the p\ sequence models, the number of superpatterns with 
unique expressions will not grow after considering order higher than Oj. The number of 
the superpatterns with unique expressions for all p\ sequence models will not grow after we 
consider the order higher than the maximum value of Oj , for j = 1 , . . . , p\ . If this maximum 
value is smaller than p, the number of the compressed parameters converges before considering 
the highest possible order, p. On the contrary, the number of the original parameters (the 
regression coefficients for those interaction patterns expressed by some training case) will 
keep growing until considering the highest order, p. 




that is also expressed by x*. If A t ^ x* t for any t with I t = 1 and A t ^ 0, 



in (3.40). 
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3.5.2 Experiments Demonstrating Parameter Reduction 

In this section, we use simulated data sets to illustrate our compression method. We apply 
our compression method to many training data sets with different properties to demonstrate 
how the rate of parameters reduction and the number of compressed parameters depend on 
the properties of the data, but without running MCMC to train the models and assessing 
the predictive performance with test cases. 

We generated data sets with varying dimension p, and number of possibilities of each 
feature, K, the same for all p features. We consider varying order, O. In all datasets, the 
values of features are drawn uniformly from the set {1, . . . , K}, with the number of training 
cases N = 200. We did three experiments, in each of which two of the three values p, K and 
O are fixed, with the remaining one varied to see how the performance of the compression 
method changes. The values of p, K and O and the results are shown in Figure 3.17. 

From Figure 3.17, we can informally assess the performance of our compression method 
in different situations. First, when p and O are fixed, as shown by the top plots, the num- 
ber of compressed parameters decreases with increasing K, but the number of the original 
parameters does not, showing that our compression method is more useful when K is large, 
in other words, when K is larger, more patterns that do not need to be represented explic- 
itly will exist. Note, however, that this does not mean the predictive performance for large 
K is better. On the contrary, when K is larger, each pattern will be expressed by fewer 
training cases, possibly leading to worse predictive performance. Second, as shown by the 
middle plots where O and K are fixed, the numbers of both the original parameters and 
the compressed parameters increase very quickly with p, but their ratio decreases with p. 
In the bottom plots with fixed p and K, the number of the original parameters increases 
with O but at a much slower rate than with p. The number of compressed parameters have 
converged when = 4, earlier than = 7. 
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Figure 3.17: Plots of the number of the compressed parameters (the lines with o) and the 
original parameters (the lines with x ), in log scale, their ratios (the lines with A) for Bayesian 
logistic classification models. The number of training cases is 200 for all data sets. The titles 
and the horizontal axis labels indicates the values of p, K and O in compressing parameters, 
where p is the number of features, K is the number of possibilities of each feature, and O is 
the order of interactions considered . 
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3.5.3 Experiments with Data from Cauchy Models 

We tested the Bayesian logistic classification models on a data set generated using the true 
model defined in Section 3.2.3. The number of features, p, is 7, and each feature was drawn 
uniformly from {1, 2, 3}. For generating the responses, we consider only the interactions from 
order to order 4, we let a — l/o, for o — 1, . . . , 4, then generated regression coefficients /3's 
from Cauchy distributions, as shown by (3.3), except fixing the intercept at 0. We generated 
5500 cases, of which 500 were used as training cases and the remainder as test cases. We 
did experiments with and without the parameters compressed, for order O = 1, . . . ,7. 

From these experiments, we see that our compression method can reduce greatly the 
number of parameters and therefore saves a huge amount of time for training the models with 
MCMC. The number of the compressed parameters does not grow any more after considering 
= 4, earlier than = 7. After compressing the parameters, the quality of Markov chain 
sampling is improved, seen from Figure 3.19, where 5 out of 6 experiments show that the 
autocorrelation of the a Q decreases more rapidly with lag after compressing parameters. The 
predictive performance with and without the parameters compressed are similar, with the 
optimal performance obtained from the Cauchy model with order O = 4, as is expected 
since this is the true model generating the data set. From the left plot of Figure 3.21, we 
see that a Cauchy model allows some (3 to be much larger than others, whereas a Gaussian 
model keeps all of j3 in a small region. For those truly small (3, Cauchy priors can keep them 
smaller than Gaussian priors can, as shown by the right plot of Figure 3.21. 

3.6 Conclusion and Discussion 

In this chapter, we have proposed a method to effectively reduce the number of parameters 
of Bayesian regression and classification models with high-order interactions, using a com- 
pressed parameter to represent the sum of all the regression coefficients for the predictor 
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Figure 3.18: Plots showing the reductions of the number of parameters and the training time 
with our compression method using the experiments on a data from a Cauchy model. The 
upper-left plot shows the numbers of the compressed and the original parameters based on 
500 training sequences for O — 1,2, ... , 7, their ratios are shown in the upper-right plot. In 
the above plots, the lines with o are for the methods with parameters compressed, the lines 
with x are for the methods without parameters compressed, the dashed lines are for the 
methods with Gaussian priors, and the dotted lines are for the methods with Cauchy priors. 
The lower-left plot shows the training times for the methods with and without parameters 
compressed. The lower-right plot shows the prediction time only for the methods with 
parameters compressed. 
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Figure 3.19: The autocorrelation plots of cr 's for the experiments on a data from a Cauchy 
model, when the order = 7. We show the autocorrelations of a , for o = 5,6, 7. In the 
above plots, "Gaussian" in the titles indicates the methods with Gaussian priors, "Cauchy" 
indicates with Cauchy priors, "comp" indicates with parameters compressed, "no comp" 
indicates without compressing parameters. 
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Error Rates with Cauchy Priors 
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Figure 3.20: Plots showing the predictive performance using the experiments on data from a 
Cauchy model. The left plots show the error rates and the right plots show the average minus 
log probabilities of the true responses in the test cases. The upper plots show the results 
when using the Cauchy priors and the lower plots shows the results when using the Gaussian 
priors. In all plots, the lines with o are for the methods that compress parameters, the lines 
with x are for the methods do not compress parameters. The number of the training and 
test cases are respectively 500 and 5000. The number of classes of the response is 2. 
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Comparison of Medians of s Comparison of Medians of s (Larger Scale) 
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Figure 3.21: Scatterplots of medians of all f3 of the last 1250 iterations of Markov chain 
samples, for the models with Cauchy and Gaussian priors, from the experiment with data 
from a Cauchy model, with the order = 4, and with the parameters compressed. The 
right plot shows in a larger scale the rectangle (—1,1) x (—1,1). 

variables that have the same values for all the training cases. Working with these com- 
pressed parameters, we greatly shorten the training time with MCMC. These compressed 
parameters can later be split into the original parameters efficiently We have demonstrated, 
theoretically and empirically, that given a data set with fixed number of cases, the num- 
ber of compressed parameters will have converged before considering the highest possible 
order. Applying Bayesian methods to regression and classification models with high-order 
interactions therefore become much easier after compressing the parameters, as shown by 
our experiments with simulated and real data. The predictive performance will be improved 
by considering high-order interactions if some useful high-order interactions do exist in the 
data. 

We have devised schemes for compressing parameters of Bayesian logistic sequence pre- 
diction models and Bayesian logistic classification models, as described in Section 3.4 and 3.5. 
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The algorithm for sequence prediction models is efficient. The resulting groups of interaction 
patterns have unique expressions. In contrast, the algorithm for classification models works 
well for problems of moderate size, but will be slow for problems with a large number of cases 
and with very high order. The resulting groups of interaction patterns may not have unique 
expressions, requiring extra work to merge the groups with the same expression afterward. A 
better algorithm that can compress the parameters in shorter time and have a more concise 
representation of the group of parameters may be found for classification models, though the 
improvement will not shorten the training time. 

We have also empirically demonstrated that Cauchy distributions with location parameter 
0, which have heavy two-sided tails, are more appropriate than Gaussian distributions in 
capturing the prior belief that most of the parameters in a large group are very close to 
but a few of them may be much larger in absolute value, as we may often think appropriate 
for the regression coefficients of a high order. 

We have implemented the compression method only for classification models in which the 
response and the features are both discrete. Without any difficulty, the compression method 
can be used in regression models in which the response is continuous but the features are 
discrete, for which we need only use another distribution to model the continuous response 
variable, for example, a Gaussian distribution. Unless one converts the continuous features 
into discrete values, it is not clear how to apply the method described in this thesis to 
continuous features. However it seems possible to apply the more general idea that we need 
to work only with those parameters that matter in likelihood function when training models 
with MCMC, probably by transforming the original parameters. 
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